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Abstract 


The ability to compute rarefied, ionized hypersonic flows is becoming more important as 
missions such as Earth reentry, landing high mass payloads on Mars, and the exploration 
of the outer planets and their satellites are being considered. Recently introduced 
molecular-level chemistry models that predict equilibrium and nonequilibrium reaction 
rates using only kinetic theory and fundamental molecular properties are extended in the 
current work to include electronic energy level transitions and reactions involving 
charged particles. These extensions are shown to agree favorably with reported transition 
and reaction rates from the literature for near-equilibrium conditions. Also, the 
extensions are applied to the second flight of the Project FIRE flight experiment at 1634 
seconds with a Knudsen number of 0.001 at an altitude of 76.4 km. In order to 
accomplish this, NASA’s direct simulation Monte Carlo code DAC was rewritten to 
include the ability to simulate charge-neutral ionized flows, take advantage of the 
recently introduced chemistry model, and to include the extensions presented in this 
work. The 1634 second data point was chosen for comparisons to be made in order to 
include a CFD solution. The Knudsen number at this point in time is such that the 
DSMC simulations are still tractable and the CFD computations are at the edge of what is 
considered valid because, although near-transitional, the flow is still considered to be 
continuum. It is shown that the inclusion of electronic energy levels in the DSMC 
simulation is necessary for flows of this nature and is required for comparison to the CFD 
solution. The flow field solutions are also post-processed by the nonequilibrium 
radiation code HARA to compute the radiative portion of the heating and is then 
compared to the total heating measured in flight. 
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Section 1 
Introduction 


The aerothermodynamic environments experienced by hypersonic vehicles are characterized 
by strong shock waves, which produce high temperatures resulting in high heating of the 
vehicle. The severity of the environment is dependent upon the mission profile and the 
vehicle configuration[l]. A high-speed planetary-entry blunt-body vehicle, such as the 
Galileo probe entering Jupiter or a high speed Earth return capsule, encounters an 
environment dominated by high enthalpy effects, which include elevated levels of ionization 
and may result in significant radiative heating. The aerodynamic and aerothermal loads on a 
hypersonic vehicle set the performance requirements for various sub-systems, such as the 
thennal protection system (TPS) and the control system. The design and safe operation of 
these vehicles, therefore, require adequate definition of the flight environment. 

During the design phase, the aerothermodynamic perfonnance of a hypersonic vehicle is 
defined by use of either experimental or computational tools. However, as in most flight 
regimes, these definitions cannot be solely obtained from ground test facilities, as no facility 
can reproduce all aspects of the flight environment. This limitation is particularly true for the 
hypersonic flight regime because of the high energy requirement needed to create a true 
hypersonic environment at a reasonable scale on the ground, and especially applicable to the 
transitional/rarefied portion of the trajectory because of the very low density requirement. 
Therefore, designers rely on computational tools for the prediction of the aerothermal loads 
experienced by vehicles entering or re-entering the atmosphere from space. 

In rarefied flows, because of the breakdown of the continuum assumption, kinetic methods 
are employed. Since its introduction in the early 1960s, the direct simulation Monte Carlo 
(DSMC) method of Bird[2] has become the de facto standard for simulating low-density gas 
dynamics. For these types of flows, traditional computational fluid dynamics (CFD) methods 
are invalid because assumptions made in developing the differential equations on which they 
are based break down under rarefied conditions. In contrast to continuum CFD methods, the 
DSMC method performs a direct physical simulation of the gas at the molecular level that 
closely mimics the Boltzmann equations. In the simulation, particles are tracked in space and 
time, accounting for both gas surface interactions and intermolecular collisions. 

In the upcoming years, the ability to compute rarefied, ionized hypersonic flows will become 
increasingly important. Missions such as Earth (re)entry, landing high mass payloads on 
Mars, and the exploration of the outer planets and their satellites will continue to push the 
state of the art in DSMC simulations. These missions may result in shock-layer temperatures 
reaching into the tens-of-thousands of degrees, where the treatment of electronic energy 
levels and ionization processes becomes increasingly important. The proper treatment of 
these mechanisms is needed to accurately predict overall shock wave thickness and structure, 
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as well as heat transfer to the vehicle. Even when the reactants are minor species, the overall 
flow field is not affected by details of the energy inodes and reaction model, but a particular 
radiation signature may be strongly affected [3]. Either way, the accurate treatment and 
prediction of electronic energy distribution and ionization is an important component of the 
design of vehicle sub-systems. 

1.1 Objectives of this Work 

Although there has been much recent interest in the modeling of ionization processes [4-9], 
the inclusion of electronic energy levels in DSMC simulations has largely been overlooked 
and they are not included in many DSMC implementations. However, recent advances in 
DSMC algorithms [10- 14], specifically the Quantum-Kinetic (Q-K) Chemistry Model, have 
allowed the prediction of equilibrium and non-equilibrium reaction rates using only kinetic 
theory and fundamental molecular properties (i.e., no macroscopic reaction-rate infonnation). 
Models such as these are important especially when reaction rate data are limited or the gas is 
far from equilibrium and an engineering approximation is required. Thus far, these new 
models have only covered neutral-species reactions. It is the purpose of this thesis to expand 
the Q-K model to include electronic energy level transitions and also to include reactions 
involving charged species. 

In addition to ionization, radiation can become an important factor depending on the mission 
and thermal protection system. For example, inflatable ballutes can be used to decelerate 
spacecraft while still high in the atmosphere. The material used for the skin of the ballute 
would have a very low tolerance for heating, necessitating the knowledge of what the 
radiative heating load would be. Several methods of calculating radiative heating for 
rarefield flows[15-18] have been developed. For the current study, however, an uncoupled 
methodology used in many computational fluid dynamics radiation evaluations [19, 20] was 
adopted. While this study does not improve on the methodologies used for calculating the 
radiative heat transfer, having separate species non-equilibrium electronic temperatures 
would make radiative predictions much more accurate. 

The models developed herein will be validated against reaction rates reported in the literature 
and applied to the second flight of the Flight Investigation Reentry Environment (FIRE II) 
Project. This was an Apollo-era flight experiment to measure the radiative and convective 
heating on spacecraft material during atmospheric reentry at lunar return speeds. The 
trajectory point at 76.42 km was chosen for the computations because, while the flow is 
nearly in the transitional regime and DSMC is more appropriate, computational fluid 
dynamics (CFD) has frequently been applied at this altitude with success. Therefore, the 
DSMC simulations presented in this thesis will be compared to CFD solutions and flight data 
from this trajectory point. However, while the methods developed herein are completely 
applicable to higher altitudes, the use of CFD is not recommended because of the break down 
of the assumption of continuum flow. 

1 .2 Outline of Thesis 

The remainder of this thesis is organized as follows. The second Section provides an 
overview of the governing equations that are the basis for the simulation of fluid flows from 
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the continuum through rarefied regimes. Also, a description of the methodology and 
structure of the DSMC algorithm, both “traditional” DSMC and the more recent algorithmic 
advances, is presented as well as previously existing electronic energy level treatments and 
ionized gas models. Section Three describes the extensions to the Q-K model and their 
implementation. In Section Four, these extensions are applied to an adiabatic heat bath and 
sampled transition/reaction rates are compared to those rates found in the literature. In 
Section Five, these extensions are then compared against previous chemical reaction models 
as applied to the second Flight Investigation of Reentry Environment (FIRE II) 
experimental]. Finally, Section Six summarizes the results of the thesis, the contributions 
of the thesis to the state-of-the-art, and recommended future research directions. 
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Section 2 

Simulation of Rarefied and Weakly Ionized Reentry Flows 


In this Section, the theoretical basis of particle simulation is first summarized. Next, a 
description of the direct simulation Monte Carlo numerical method is given, including the 
traditional methods and more recent algorithmic advances, as well as previous attempts to 
include electronic energy levels and models of ionized flow. Throughout this description, 
mention will be given of any modifications that were required to be perfonned on the DSMC 
algorithm of choice. 

2.1 Theoretical Basis of Particle Simulation 

The particle simulation method has been developed independently by plasma physicists and 
aerothennodynamisists. Plasma physicists, who simulate charged particles, have developed 
the particle-in-cell (PIC) method based on the theory of charge assignment and force 
interpolation. This theory has made it possible to obtain smooth distributions of charge 
density and current density in Maxwell’s equation of electromagnetic fields. The details of 
the PIC method can be found in the classical works of Birdsall and Langdon[22] and 
Hockney and Eastwood[23]. On the other hand, aerothermodynamicists have been 
developing a particle simulation method of neutral species beginning with the work of 
Bird[2, 24], This method is called the direct simulation Monte Carlo (DSMC) method. 

Particle simulation of both charged and neutral species is necessary to examine the structures 
of high-speed, rarefied, reacting flows. However, for flows where the plasma density is 
considered low (-10%), the assumption of charge neutrality is generally assumed and the 
charged particles are treated as neutral particles[25-30]. More will be discussed on these 
techniques in Section 2.2.4. The current section describes how the probability laws used in 
particle simulation may be derived[3 1] from a relevant Boltzmann equation for a weakly 
ionized plasma and simplifications will be pointed out for our assumption of charge 
neutrality. 

2.1.1 The Boltzmann Equation 

The theoretical basis of particle simulation can be explained using a weakly ionized plasma 
consisting of electrons (e), ions (A), and molecules ( B ). If Coulomb collisions e-e, e-A, and 
A- A can be disregarded, only the remaining e-B, A-B, and B-B collisions need to be 
considered. The particle motions of species e, A, and B are governed by the respective 
Boltzmann equation for each of these species. Except for Coulomb collisions, the Boltzmann 
equation for electrons or ions is linear, but the Boltzmann equation for molecules is nonlinear 
because of the B-B collisions. Although the PIC method was developed for charged species 
and the DSMC method for neutral species, the idea behind the two methods is similar. The 
Boltzmann equation of ions, A, is written as: 

^ + v j^ + ^^(FjfA) = i\[fA{?)f(&)-fA{v)f{w)]8<r AB dndw ( 2 . 1 . 1 ) 

j A j 
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Here, Coulomb collisions are neglected and only A -B collisions are considered, E } is the force 

acting on particle A with a mass of mA,f. 4 (v ) and fs{ w ) are the velocity distribution functions 
of species A and B , v' and w' are the post-collision velocities, g - |v — vt>| is the relative 

speed, <y 4B is the differential cross section, and dll (=sin^ d% d y/) is the solid angle. In the 
remainder of the derivation, f A (v,x,t ) and f B (v,x,t) are expressed as / A (v) and f B (v ) in 

the right hand side of the equation for simplicity. The differential cross section <J 4B is a 
function of g and the deflection angle %. The force F is given by: 

F = q A (E + vxB) (2.1.2) 

where q A is the charge of the particle A and E and B are the electric and magnetic fields, 
respectively. However, for our purposes, the force experienced by the charged particles is 
effectively zero since there is an assumption of no charge separation (i.e., charge neutrality). 
For the sake of completeness, though, this tenn will be kept in the remainder of the 
derivation. 

If the method to obtain f A (v,x,t) is given, the solution of /a can be found at any time by 

step-by-step calculations started from an initial velocity distribution function. For a small At, 
we have: 


f A ( v,x,t + At) = f A (v,x,t) + At 


V dt j 


Substituting df A / dt of Equation (2.1.1) into Equation (2. 1 .3) gives us: 
f A ( v,x,t + At) = (l + At/)(l + A tD)f A (v,x,t) 
where D and J are operators defined by 

df, i a 


£>(/,)= 


-V, 


dxj m A dvj 


( F jf A ) 


J (/a ) = JJ [/a {v')fs i w ') - f A ( v)f B (w)] g(J AB dEldw 


(2.1.3) 

(2.1.4) 

(2.1.5) 

(2.1.6) 


D and / denote non-collisional and collisional operators, respectively, and therefore, 

Equation (2. 1 .4) decouples the non-collisional and collisional parts of the motion. Equation 
(2. 1 .4) is equal to Equation (2. 1 .3) if the term of order (At) 2 is disregarded such that the 
solution obtained from Equation (2. 1 .4) is first-order accurate. In general, At should be much 
less than the mean free time x (the mean free path divided by the mean speed of the 
particles). This condition should be strictly satisfied in the case when particle A is charged 
[32]. Note, however, that if particle A is a molecule, the condition can be replaced by a 
weaker one of At < x. The reason why the weaker condition can be applied is that molecular 
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trajectories are straight and, therefore, accurate for any At and that the state of the gas 
changes little in the mean free time. 


Equation (2.1.4) suggests that a method of solving the Boltzmann equation is to first obtain 
f* A = (1 + AtD) f A (non-collisional) and then obtain (l + A tj) f A (collisional). The first 
procedure is equivalent to solving: 


dt 


+ V; 



i a 


m A dv 


( F +)=° 


(2.1.7) 


by use of the initial distribution function f A (v,x,t ) . If the solution at t+At is f A (v,x,t ) , then 
the solution is f A (v,x,t) = f A [v -(P l m A ) At,x - vAt,t ) , and therefore 


/I (v + Av,x + A x,t) = f A (v,x,t) 


(2.1.8) 


where A v = i y F I m A ) At and At = vAt . Equation (2. 1 .8) implies that the particles of species 
A move according to the equation of motion 


CIV / - _ -\ 

™ A — = q A (E + vxB) 

cl [ (2.1.9) 

dx 

— = v 
dt 


These equations of motion are used in PIC simulations which means that the Boltzmann 
equation without the collision term is that solved by the PIC method. When Equations 
(2. 1 .9) are solved for all particles, which is trivial for DSMC since there is no net force, we 

have f A (v,xj). The second procedure is equivalent to solving 


= JJ[/a ( v ')fs ( w ' ) - Ia ( v )f B ( w)]go AB dQdw 


(2.1.10) 


using the initial distribution function f A (v,x,t ) . The solution of Equation (2. 1 . 10) at t+At is 
f A (v,x,t + At ) . 

Equation (2.1.10), which does not include the term df A / r).v ; , represents the collisional 

relaxation of the velocity distribution in a spatially uniform state of species A. Equation 
(2.1.10) can be applied only to a small cell where the gradient of density, temperature, and 
flow velocity of species A can be regarded as uniform. Because of this requirement, the cell 
dimension should be of the order of the mean free path of species A, and is the basis for the 
cell size of the DSMC method. 

A solution of equation (2.1.10) can be obtained by an equivalent stochastic process[33]. 

First, in particle simulation, it is sufficient to consider the motions of a set of particles 
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randomly sampled from all particles. The sampled particles are referred to as simulated 
particles. The number of simulated particles in a cell is much smaller than is the number of 
real particles. The weight, W, is defined as the number of real particles represented by each 
simulated particle. After moving all simulated particles by At, we identify the simulated 
particles in each cell by examining the positions of the particles. Let N A be the number of 
simulated particles in a cell and v A] ,v A2 . . be the particle velocities. The velocity 

distribution function f A (v,x,t) for a cell located at x is expressed as 

f A {v,x,t) = ^-^5 i (v-v A ) (2.1.11) 

where n A is the number density of the cell, S 3 ( ) is the Dirac delta function, and v Al is the 
solution of Eq. (2. 1 .9) at time t+At. The function <5 3 (v - v Ai ) is the probability density 
function for the velocity v Aj of particle A h Similarly, the velocity distribution function of 
species B at the end of the particle’s motion is 

f B {w,x,t)=^^S 3 (w-v B ) (2.1.12) 

yV B j = 1 

where v Bi is the velocity of simulated particle B,-. The solution of f A (v,x,t + At) can be 
obtained from 

f A (v,x,t + At) = (1 + A tj)f* (v,x,t) 

= fl (v,x,t) + At JJ[/ A (v)f B (w) - f A (v)f B (vv)]ga'' B dQdw 

Substitution of Equations (2.1.1 1) and (2. 1 .12) into Equation (2. 1.13) yields 

f A (v,x,t + At) = -^-^ j f Ai {v) 

^ A i= 1 

where f Ai (v) is the probability density function for the velocity of particle A t at t+At. It is 
given by 


(2.1.13) 


(2.1.14) 


f A , (v ) = (l - P Ai )<$ 3 (v - v A ) + P A Q Ai (v ) (2.1.15) 

The structure of Eq. (2.1.15) states that P A i is the collisional probability of particle^/. Note 
that if P A i = 0, f M (v) = 8 3 (v - v A; .) , and there is no change in velocity because of a collision. 
The probability P Ai is given by 




j = i 


(2.1.16) 
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where 


^4, , Bj N B n B8A i ,B j ®T (<? Ar B ; )^ 


(2.1.17) 


Here, g Ai Bj = |v A; . - v Hj is the relative speed and (j ' 11 is the total collision cross section. 

Equation (2. 1.16) reveals that P A i,Bj is the probability of particle A, colliding with particle B j 
in At. The function Q Aj (v) is the probability density function for the post-collision velocity 
v . It is written as: 


Qa , » = 



(2.1.18) 


Here, the ratio P A i,Bj/PAi represents the conditional probability that the collision partner of 
particle A,- is particle Bj under the condition that particle A, collides in At. Here, Q Ai ,bi (^) is 

the probability density function for the post-collision velocity of particle A, when its collision 
partner is particle Bj, which is given by: 


where 


Qa„Bj{ v ) 


g AB (8a,j 1j ’X)8§8'\-8 Ai j 1j ) 

(g A ., B ,) 2 M 3 B o* B (g Ai , B J 


|' = — (v- 


M AV A - M BV Bi 


(2.1.19) 


(2.1.20) 


M a = m A /(m A +mB), M B = mB/(m A +niB), X is the angle between g' and v Ai -v Bj , and 
& 4B (gAi,Bj,%) is the differential cross section. From Eq. (2.1.19), we have 


Q A .B i ( v ) dv = s ( Isl 



0 AR {g Ai . Hr x)s™xdxdV 

a T B (gA„B,) 


(2.1.21) 


where dv = M 3 B dg' = M 3 B g' 2 dg sin %d%dii/ . Note that M 3 B is the Jacobian of the 
transformation from v to g ' , and then Eq. (2. 1 .20) yields the post-collision velocity v Ai (= v) 
of particle A t 

v' Al = 1 ~ (m A v A + m B v B + m B g') (2. 1 .22) 

m A + m B v 1 ’ 
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where g'(= v' Ai - v' Bj ) is the post-collision relative velocity and Ig'l = g ALBj . The probability 
that g' is included in the solid angle sin jr/jr/y/ is <J AB (gAi'Bj’X^XdXdV / erf (g M<Bj ) . 
The direction {XW) of 8 can be sampled using this probability. 

Equations (2.1.16) and (2.1.17), which are derived from the Boltzmann equation, can also be 
obtained by use of the elementary free-path theory. The collision frequency of particle Aj is 


v 


A, , 



where g. 4 i,B is the relative speed between^, and some particle B, and the bar denotes the 
average of particles B. Therefore, the average of all simulated particles B in a cell and the 
collision probability of particle Aj take the form 


= V A, At = 


n B At 

^7 



j = i 



This equation agrees with Equation (2. 1.16). We cannot obtain the probability laws to 
detennine the post-collision velocities from the free-path theory, however. In many cases, 
we can write such probability laws intuitively because the laws derived from the Boltzmann 
equation agree with the physical image of the collision. 

The probability laws used to detennine the velocity of particle Aj at t+At are then: 

1) Particle Aj collides in At with a probability of Pau 

2) Under the condition that particle Aj collides, its collision partner is particle Bj with a 
probability of P.u.m Pu- 

3) If the collision partner is Bj, the post-collision velocity v' Ai is given by Equation 

( 2 . 1 . 22 ). 

Therefore, the particle simulation method can be derived from the Boltzmann equations of 
ions and is true for any species. Although the Boltzmann equation for a simple gas is 
nonlinear because of collisions between like molecules, the simulation method was derived 
in 1980[34], Before this work, those in the field of rarefied gas dynamics had regarded the 
particle simulation method as a kind of numerical experiment. By use of the measure theory, 
mathematicians[35] verified that Nanbu’s particle simulation method is an appropriate 
method for solving the Boltzmann equation. 

The most important ideas of the particle simulation method can be summarized as: 1) 
decoupling of collisionless motion and collisions; and 2) calculating collisions independently 
in each cell. These basic concepts have been shown to naturally result from the operator 
splitting in the fonn of ( \+AtJ)(\+AtD ). Now consider the conditions for choosing At. The 
operator splitting is allowed if and only if At|/(/ A )| «: f A and At|D(/ A )| <sc f A . The first 
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condition states that the effect of collisions in At must be small; therefore the collision 
probability P Ai — At / r is small, where x is the mean free time. The time step must satisfy 
the condition At <SC T . The second condition requires that | Ax <sc |.r and |Av| <sx |v| in Eq. 
(2.1.8). That is, VAt <sc x and (F / m A ) At v , where v is the mean particle speed and 
F / m A is the mean acceleration. The condition that v At x is usually replaced by the 
Courant conditionvAt < A c , where A c denotes the cell size. Therefore, a particle’s 
displacement in At must be less than the cell size. The condition (f / m A ) At <sc v requires 

that the change in the particle’s speed in At must be much less than the original speed for PIC 
simulations. Therefore, a smaller At should be used for an external field. 

2.2 Traditional Direct Simulation Monte Carlo Methods 

The DSMC technique was developed primarily by Graeme Bird, now Emeritus Professor at 
the University of Sydney, Australia, nearly a half century ago. Since its inception, this 
technique has been used in many different applications ranging from the simulation of micro- 
electromechanical devices (MEMs), to various industrial and phannacological deposition 
simulations, to the simulation of (re)entry and other rarefied flow fields. 

In the work outlined herein, the DSMC method is used to simulate the Boltzmann equation 
for rarefied/transitional reentry flows in which the Navier-Stokes equations are not valid 
throughout the entire flow field. Specifically, the algorithm employed is the DSMC Analysis 
Code (DAC)[36, 37], which was developed by a team led by LeBeau beginning at the NASA 
Langley Research Center and continuing at the NASA Johnson Space Center. The DSMC 
algorithms in the release version of DAC are based on those developed by Bird[2], 

Regardless of the DSMC code used, the basic structure of the algorithm is the same. The 
following is included to briefly describe the flow of a generic DSMC algorithm. 

2.2.1 Collision Algorithm 

In each computational cell, particles are chosen to collide in such a way as to reproduce the 
required macroscopic collision rate. The simplest implementation of collision pair selection 
is to disregard particle position within the cell when selecting the pair. This methodology 
restricts the size of a computational cell in that it must be smaller than the local average mean 
free path of the flow field. This was the first selection methodology and is still utilized 
today. The next selection methodology was to select the first particle randomly in the cell 
and then select the particle closest to the first excluding its previous collision partner, or the 
so-called virtual sub-cell (VSC) algorithm of LeBeau[38], This methodology allows for a 
larger collision cell size as long as the resulting mean collision distance divided by mean free 
path remains at least as small as the resulting value from the original selection methodology. 
Other selection methodologies have been proposed[39, 40], but the two discussed above are 
the most prevalent. The VSC algorithm is that which was used in the current study. 
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In many modern DSMC codes, Bird’s No Time Counter (NTC) method is used to select 
colliding pairs[2]. At each time step and in each computational cell, the number of particle 
pairs to be tested for collisions is calculated as follows: 

N, M =\nN(o T g)^to (2.2.1) 


where n is the number density of gas particles in the cell, N is the average number of 
simulator particles in the cell, At is the simulation time step, g is the relative velocity of a 
colliding pair of particles, and Or is the total collision cross section of the colliding pair. 

Each of these pairs are then selected to collide if the probability: 

p = (2 2 2) 

coll / \ 

( o T g) 

\ i ° / max 


is less than a random number uniformly distributed between zero and one. This method 
assumes that the collision rates between species will be implicitly reproduced, even though 
all species are lumped into the same group for the pair selection process. This method can 
become inefficient if the product of the collision cross section and relative velocity becomes 
very much different between species pairs, as is the case when electrons are introduced into 
the flow, because the maximum value of o T g is very large and results in many null 
collisions. This is the basic algorithm used by the DAC code. 

In order to address this issue, the collision algorithm in DAC was modified as follows[31]. 
The collision algorithm can be divided into two groups. One includes collisions between 
similar particles, and the other collisions between dissimilar particles. The collision tenn of 
the Boltzmann equation is nonlinear for like collisions and linear for unlike collisions. 
Consequently, the collision algorithm for like collisions is different from that for the unlike 
collisions. 


There are five classes of collisions involving charged particles: electrons with neutrals, ions 
with neutrals, electrons with ions, ions with ions, and electrons with electrons. The 
interactions of the last three classes obey the Coulomb force law: 


p _ 1 

4tt£ 0 r 2 


(2.2.3) 


where e 0 is the permittivity of free space and r is the distance between the particles. 

Bird[41] and Taylor el. al. [25], however, used the variable hard sphere (VHS) model of Bird 
to describe the first four classes of collisions by setting the diameter of the heavy ions equal 
to the diameter of their neutral counterparts. The reference diameter of the electrons is set 
equal to le-10 m, which is orders of magnitude larger than the classical value of 6e-15 m. 
This choice was somewhat arbitrary, but Bird did mention that the electron diameter was 
varied in his simulations by a factor of three with no significant effect on the flow field 
parameters. Neither modeled electron-electron collisions. This is the approach that will be 
followed herein. 
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In general, then, one can divide the collision algorithm into two for weakly ionized plasmas, 
neglecting Coulomb collisions: 1) collisions between unlike particles; and 2) collisions 
between like particles. The algorithm for type (1) is derived from the linear Boltzmann 
equation as in Section 2.1.1, and that for type (2) is from the nonlinear Boltzmann equation. 

Let us consider the elastic collisions: 

A(v a ) + 5 (v b )^A(v;) + 5(v;) 

If gas B is at rest and in equilibrium at temperature T B , the probability that particle A 
elastically collides with some particle B is: 


P cou = n B8° T {g) At (2.2.4) 


where n B is the number density of gas B, g = | v A - v B is the relative velocity, and the average 
is taken over the Maxwellian distribution of v B . The average in Equation(2.2.4) is for all 
simulated particles in the cell, which gives: 


P coU=^'Zgj°T{gj) A < 


(2.2.5) 


B 7=1 


where Nb is the number of simulated B particles. If the collision is treated by the hard-sphere 
model, the scattering is isotropic. The post-collision velocities are then: 


v\ 


ttt( 


m A + m B 


m A V A + ™B V B + m B V A 


Vo = 


— ( 

W I ' 


m A + m B 


m A V A + m B V B ~ m A \ V A ~ V B 


~Vb\R) 

*) 


(2.2.6) 


where the two R's are the same unit vector with a random direction. However, the hard 
sphere model is generally insufficient for most applications because it does not match the 
gas’s viscosity or diffusion rate as a function of temperature. On the other hand, the most 
accurate Lennard-Jones potential is not suited for the DSMC method because the collision 

calculation is very time-consuming. But, the simple inverse-power potential V(r)= a/r a is 
sufficiently accurate and very efficient for collision calculations. The differential cross 
section a for the inverse-power potential takes the form: 

cr = Ag~ 4/a (2.2.7) 


where g is the relative speed and A is a function of the deflection angle X- If A is independent 
of x, the probability a dQ./ a T that the post-collision relative velocity is directed into the 

solid angle dQ becomes dQ/4n, or the scattering is isotropic and the post-collision velocities 
are the same as that of the hard-sphere model. The variable hard-sphere (VHS) model was 
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proposed by Bird[41] to avoid the divergence of the total collision cross section and using 
isotropic scattering for the inverse-power potential. In the VHS model, the coefficient A is 
constant and, once a gas is selected, A and a can be chosen to best tit the Chapman-Enskog 
viscosity to experimental data. This leads to a power law temperature dependence of the 
coefficient of viscosity such that: 


ii r r0) 

/i 1 


where 


(O = 


a + 4 
2a 


(2.2.8) 


(2.2.9) 


and to is bounded by 0.5 and 1.0 at the hard-sphere and Maxwell molecule limits, 
respectively. The VHS parameters used in the current study are discussed in Section 5.3. 

2. 2. 1.1 Simple Gas 

For a single species gas in a computational cell, let v x ,v 2 ,...,v N be the velocities of the 

simulated particles. The probability P t that molecule i collides in At with some other particle 
is given by Equation (2.2.5). For a simple gas, this is: 

F i =j;lg,MsM , =t p » < 2 - 2 - 10 > 

n j = i j = i 


where g tj = v ( . - Vj and 


p ij = N~'ng ij o T (g 0 . ) A t (2.2. 1 1) 

Note that P u = 0 because gu = 0. The tenn for j = i should therefore be omitted in Equation 

( 2 . 2 . 10 ). 

If gy is then replaced by g max , the resulting probability is: 

(^L (2.2.12) 

This probability is greater than Pij. For Pj = ( Pij ) max , Equation (2.2. 10) takes the form: 

(2.2.13) 

The number of collisions in At for ( P,) ma x is 

^max = iZWmax = ^{N ~ l) g Jmx (J T (g^) At (2.2.14) 


13 



The factor of V 2 on the right hand side of Equation (2.2.14) is to avoid double counting 
collisions between like particles. N max is called the maximum collision number. The 
computing time necessary to obtain the real collision number 

N„,, = 1'L'Lp„ (2.2.15) 

1=1 j = 1 

is proportional to N 2 . This is not insignificant because collisions in all cells in a 
computational domain must be calculated at each time step. In the null-collision method, the 
task is avoided by regarding N max as the collision number in At. To compensate for excess 
collisions, N ma x-N c , some of the N max collisions are regarded as null collisions. We now 
define: 



where q tj = g tj C 7 r (g,. ( j / g max <7 7 . (g max ) . This means that the pair (/,/) first makes a tentative 

collision with a probability of ( Pij) m ax ■ The tentative collision results in a real collision with a 
probability of qj and a null collision with a probability of 1 -c/ (/ . The number of tentative 
collisions is ( Pij) max times the number of pairs N(N- 1 )/2; the product is equal to N max . A 
tentative collision pair can be chosen randomly because (Pij) m ax is common to all (/,/'). So, to 
summarize, we first obtain N max . In general, N max is not an integer. Let p = iV max - [yV max J , 

where |_^ T „ax J is the integral part of N max . We then set A max = \_N m ax J + 1 with a probability 
of p and N mm = [iV max J with a probability of 1 -p. The following two steps are then repeated 
N max times: 

1 . Choose a pair (ij) of i 4- j randomly and obtain q tj . 

2. Call a random number U. If U < qij, the pair collides. If they do collide, the pre- 
collision velocities are replaced by the post-collision velocities given in Equation (2.2.6) 
where the particles A and B are of the same kind. 

2. 2. 1.2 Gas Mixture 

For the sake of simplicity, consider a binary mixture of species A and B. The types of 
collisions are A -A, B-B , and A-B. The A-B collision is a collision between unlike particles. 
The other two types are treated as in the gas of a simple gas above. Let Na and Nb be the 
numbers of simulated particles A and B in the cell. The probability that molecule A, collides 
with some molecule B in At is 


p Ai =n B8^T BAt (2.2.17) 

where g is the relative speed and <Jj B is the total cross section for the A-B collision. The 
average is for all simulated molecules B in the cell, so the probability takes the form: 
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(2.2.18) 


p = y p 

A Z-i A - B i 
j = i 


where 




(2.2.19) 


If we again introduce the idea of the maximum collision number where gAija is replaced by 
Smax •> we have the upper bound of Paw'- 

PZ=N;'n»sZ<(sZ)A l ( 2 . 2 . 20 ) 


The maximum number of A-B collisions in At for P^ x is 


n a n a n b 

p AB = y V p AB 

max g—i g—i max 

»= 1 <=1 j = 1 


=y 

max 


(2.2.21) 


Note that the factor of Vi is not necessary for A-B collisions. Substitution of Eq. (2.2.20) into 
Eq. (2.2.21) gives us: 


= n B^ASZ^T B (kmfx ) 
The number density n B is expressed as: 


(2.2.22) 


n B 


wn b 

K 


(2.2.23) 


where V c is the cell volume and W is the weight. As discussed previously, one simulated 
particle represents W real particles. Presently, assume that all the particles in a cell have a 
common weight. The collision procedure in the previous section is then repeated for all 
species combinations A- A, A-B , and B-B using the updated definitions of number of 
collisions and probabilities for like and unlike collision pairs as required. 

2.2.2 Internal Degrees-of-Freedom 

A typical application of the DSMC method involves the calculation of many billions of 
intennolecular collisions, and this leads to the requirement for computationally simple 
molecular models. This has led to the development of “phenomenological” models that 
focus on the generation of accurate results for the observed gas properties but do not attempt 
to provide the most realistic model of individual collisions at the microscopic level[13]. The 
developers of classical kinetic theory faced a similar problem in that they required 
mathematically tractable molecular models that permitted the derivation of analytical 
expressions for the macroscopic transport properties. The classical models were the obvious 
choice for early DSMC work, but were found to have serious shortcomings that have been 
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overcome by phenomenological models that have been introduced in the context of the 
DSMC method. 


We have already seen a phenomenological model be applied to elastic collisions by way of 
the variable hard sphere model, which combines the unrealistic but simple isotropic 
scattering of a hard sphere with a realistic variation of diameter with the relative speed of the 
collision pair. This model leads to accurate results for viscosity and heat conduction 
coefficients. 

In the same way, classical kinetic theory did not lead to any useful models for inelastic 
collisions. For example, the rough sphere model retains the deficiencies of the hard sphere 
model with regard to the transport properties, and it is unable to deal with the quantum 
effects that cause most gases to have fewer than three rotational degrees of freedom and has a 
fixed and unrealistically fast rotational relaxation time. This problem was solved in the 
context of the DSMC method by the introduction of the Larsen-Borgnakke, or L-B, 
model[42] for rotation which can be regarded as the archetypal phenomenological procedure. 
It is an add-on to the VHS model and, for a fraction of the collisions that are chosen to match 
measured rotational relaxation rates, the post-collision rotational energies are selected from 
the equilibrium distribution that correspond to the collision energy. The L-B model was 
extended to the vibrational modes, but it assumed that the internal energies were 
continuously distributed and was found lacking. The situation was corrected by the 
introduction[43] of the quantum version of the L-B model. 

DSMC methods test colliding particles for possible internal energy exchanges or chemical 
reactions. The collision selection rules use probabilities that must be defined such that they 
result in the desired macroscopic relaxation rate behavior. This typically involves identifying 
a macroscopic rate equation that the simulation attempts to match through appropriate 
definitions of selection probabilities P. These rate equations often contain parameters such 
as collision numbers Z, characteristic times x, or rate coefficients that are estimated from 
experiments or analysis. Ambiguity regarding specific definitions of the targeted 
macroscopic rates and coefficients has often led to confusion in defining the appropriate 
selection probabilities in DMSC methods. 

In studies of vibrational relaxation, Haas[44] developed a relationship between the 
appropriate DSMC selection probability P and the resulting Landau-Teller macroscopic 
relaxation behavior. This relationship was later generalized, allowing for rotational 
relaxation, by Lumpkin et al. [45] in the form of a “correction factor” for collision numbers Z. 
It was observed that the DSMC probabilities P, in general, cannot be specified by the 
assumption P=l/Z, which has been employed throughout the DSMC community. Dissimilar 
applications of the selection probability and incompatible definitions of collision number led 
to disagreement among DSMC researchers regarding these correction factors. 

Haas et al. [46] went on to clarify the definitions of macroscopic rate parameters employed in 
rotational relaxation, including differing interpretations of Jeans’ equation[47]. It was shown 
that relaxing energies, when plotted against the cumulative number of molecular collisions, 
lead to a single solution, regardless of the assumed intermolecular potential in an adiabatic 
bath. Differing definitions of collision time T c between experiments and DSMC were 
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clarified. The two most widely used selection methodologies were discussed, and DSMC 
rotational relaxation probabilities P, specific to each methodology, were derived to match the 
macroscopic rate equations exactly. Prevalent selection methods cannot reproduce Jeans’ 
equation exactly for the case of heteromolecular collisions, nor can they reproduce Jeans’ and 
Landau-Teller[48] rate behavior for simultaneous rotational and vibrational relaxation. Thus, 
a new selection methodology was introduced that was better suited to the general case of gas 
mixtures experiencing multimode internal energy relaxation. Although constant Z values 
were used, the developments were valid for variable Z as well. 

Gimelshein et al. [49, 50] proposed a simplified inelastic collision algorithm using particle 
selection prohibiting multiple relaxations. They felt that the selection methodology proposed 
by Haas[46] lead to an unnecessarily complicated functional dependence between the 
probability and the collision number. To avoid the complexity, a straightforward selection 
methodology was developed which also resulted in relaxation behavior matching the Landau- 
Teller equation when the correct correction factor was used. 

Gimelshein et a/.[50] took a closer look at vibrational relaxation and derived an exact 
relationship that connects the vibrational relaxation number, Z vi b° , used in the DSMC 
method and that employed in continuum simulations. An approximate expression for 
Z vi b DSMC was also derived that is cost-effective and applicable when translational temperature 
is larger than the vibrational temperature. 

2.2.2. 1 Rotational and Vibrational Energy Transfer 


Models for rotational and vibrational energy transfer with translational energy in a collision 
have been explored previously. The Larsen-Borgnakke (L-B) model[42] is usually used to 
describe the rotational and vibrational energy transfer in the DSMC method. A fraction of 
the total number of collisions, 1 /Z V7 v, or 1 /Z rot , lead to the energy exchange, where Z vi /, and Z rot 
are the vibrational and rotational collision numbers, respectively. Originally, the internal 
energy modes were treated as continuous. These energy modes are, of course, discrete. To 
address this issue, Bergemann and Boyd introduced the algorithms of the L-B model with 
discrete vibrational and rotational energy levels for diatomic molecules[43, 51]. Gimelshein 
et al. then extended the L-B model for the case of polyatomic gas mixtures with discrete 
vibrational and rotational energies[49]. 


The temperature dependent nature of rotational and vibrational relaxation rates has been 
developed by Parker[52] and Millikan and White[53], which provides the characteristic 
relaxation times as functions of the temperature for a wide range of pure and mixed diatomic 
gases. In the L-B model, the vibrational and rotational post-collision energies are modeled 
according to the local equilibrium distribution functions. The derivation of these models for 
application in the DSMC method can be viewed in Reference [2], but the vibrational 
relaxation model resulting from Millikan and White, which is used in this study, is: 


Z 


vib 



V T m J 



(2.2.24) 


where Ci and C? are constants used to fit the experimental data. 
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However, these lead to the incorrect macroscopic relaxation rate. Lumpkin and Haas et 
al. [45, 46] introduced a correction factor for the rotational relaxation rate in the DSMC 
method. Most continuum analyses and computational codes use Jeans’ equation[47] to 
model rotational relaxation. In a single species gas, Jeans’ equation uses the instantaneous 
equilibrium rotational energy E rot (/) and a time constant, which can be written as the product 
of the instantaneous collision time x c and a characteristic collision number Z as: 


dE ro, _KoX t )- E roX t ) 


(2.2.25) 


where E rot is the mean rotational energy per molecule in the ensemble, and t represents time. 
Denoting the rotational degrees of freedom by £ rot , E rot is defined by: 



kT 


(2.2.26) 


and will change with time as the translational temperature T varies during relaxation under 
adiabatic conditions. Similarly, x c may vary with T, depending upon the intennolecular 
potential assumed for the gas. Most DSMC methods employ the VHS model, which assumes 
rigid-sphere collision mechanics, but establishes the collision rate in a manner consistent 
with an inverse-power potential. The relative translational energy of colliding particles near 
equilibrium, as a result of biasing due to collision selection, is distributed with (^,-ans degrees 
of freedom given by: 


C™=4-- = 5-2 to (2.2.27) 

a 

Note that this can differ from = 3 associated with the translational energy of an 
equilibrium ensemble of particles in a reservoir. Biasing from collision selection plays a 
significant role in the relaxation of the energy in a reservoir and becomes apparent in the 
derivation of the relaxation probability. Haas[46] derived the relationship for the particle 
selection method for a single species gas (allowing for multiple relaxation) as: 


pDSMC 
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(2.2.28) 


that correctly reproduced the relaxation rate given by Jeans’ equation. He goes on to say that 
P part = 1 IZ could be a fair assumption for adiabatic relaxation, but only if one redefined Z by 
using adiabatic relaxation data that has been cast into the isothennal fonn of Jeans’ equation 
as suggested by Bird[2]. As such for adiabatic flows, P par t = ( 1 /Z ISirc /) = (l/Z)(l+£ rot /£r) can 
be a fair assumption to the exact expression. Better yet, the probabilities were shown to 
approach asymptotic limits described by: 
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(2.2.29) 


which represents an excellent approximation when Z> 10 for the particle-selection 
methodology. 

Gimelshein[50] then analyzed the vibrational relaxation rates in the DSMC method. The 
main difference between rotation-translation and vibration-translation relaxation is that while 
the rotational mode can be describes as having a constant number of degrees-of-freedom at 
all temperatures except cryogenic, the number of vibrational degrees-of-freedom varies with 
the temperature. For a simple harmonic oscillator, the number of vibrational degrees-of- 
freedom Cv.b is: 


r 

ex P (e„ t /r)-i 


(2.2.30) 


where 0 v ib is the characteristic vibrational temperature of the mode. This means that the 
relationship for rotational relaxation probability given above are not applicable to the 
vibration-translation energy transfer. Gimelshein derives the relationships that allow one to 
match the vibrational relaxation rate in DSMC with that predicted by the Landau-Teller 
equation[48]: 


dE, 


vib 


E vib {t) E vib {t) 


dt 


_ ' ycont 

^coll^vib 


(2.2.31) 


where t is again time, E vlb is the mean vibrational energy per molecule in the ensemble, and 
E vib is the instantaneous equilibrium vibrational energy, defined as: 


E-rib ~ 


C (t ) 

^ vib \ trans / 


kT 


(2.2.32) 


He goes on to derive the expression for the DSMC collision number for vibration as: 

r (t —C (t i 

vib \ trans / trans -> vib \ vib ) vib 1 


p DSMC 
vib 


Cr (E trans -r) 


(2.2.33) 


The probability of a vibrationally inelastic collision in the DSMC method is then E = 1/ Z'"'" 

. Although exact, Equation (2.2.33) is costly for direct use in DSMC computations because it 
involves the calculation of the vibrational and translational temperatures, and then the 
subsequent numerical determination of T’ (post collisional temperature). However, the exact 
expression can be significantly simplified when the vibrational and translational temperatures 
are close: 
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(2.2.34) 


The above equation is applicable only when the difference between the vibrational and 
translational temperatures is small. However, it was found that if the translational 
temperature is used as T, Equation (2.2.34) reproduced well the correct expression in 
Equation (2.2.33) for all cases where the vibrational temperature was lower than the 
translational one. 

2. 2. 2. 2 Electronic Energy Levels 

In the past, Bird[ 1 5] and Carlson et al. [17, 25] modeled electronically excited levels and 
electronic transitions in DSMC. In their work, each excited particle was assigned a 
distribution over all available states, assuming that the electronic energy mode was in 
equilibrium. Gallis and Harvey[16, 54] modeled non-equilibrium thermal radiation in 
DSMC as well as electronic excitation processes. However, for their method, accurate cross 
sections for electronic excitation and absorption were required and the uncertainty in the 
cross sections is not well known for all the species. Because of these deficiencies, a much 
more accessible non-equilibrium electronic energy model is developed in this thesis in 
Section 3.1. 


2. 2. 2. 3 Selection Methodology 

Haas[46] derives the relationship for a versatile L-B selection methodology that is applicable 
to gas mixtures and is adapted from the particle-selection method; probability, P, is applied 
to each molecule individually within the colliding pair, but the event in which both particles 
relax is prohibited (called particle-selection prohibiting double relaxation). This ensures that 
relaxation of A molecules is independent of the energy of B molecules, and the L-B 
mechanics can therefore reproduce the relaxation rate behavior for each species 
independently. The appropriate functional fonn for the relaxation probability was defined as: 
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(2.2.35) 


where the Z values are from continuum methods. This expression involves only rotational 
energy exchange. The expression for multimode energy relaxation is presented by Haas[46], 
and the reader is referred to the reference for the complete derivation. 

Gimelshein[49, 50] believed that, although correct, the above selection methodology leads to 
an unnecessarily complicated functional dependence between the probability and the 
collision number. This is especially true when multi-mode relaxation is considered. To 
avoid the complexity, a straightforward selection methodology was used that is illustrated in 
Figure 2.1 . The first step is to select a random number Rn unifonnly between 0 and 1 . 
Then, using the value of Rn, the possible relaxation events are tested, so that either rotational 
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or vibrational relaxation may occur for particle A or particle B. Multiple relaxation events 
are prohibited. 



^2 = ^r.A + Pf.B 
^3 = Pf.A + ^r,B + K* 

^r,B + Pv,A + ^v.B 


Figure 2.1: Selection methodology for inelastic and elastic collisions. 

Similar to that suggested by Haas, the simplified selection methodology allows only one 
relaxation event to take place in a collision. Vibrational and rotational relaxation rates of 
different species in the gas mixture are independent, and the necessary relaxation rates may 
also be accurately captured. It can be easily implemented and is characterized by a simpler 
dependence between the collision numbers and the relaxation probabilities. This method 
should not be applied when the total probability of internal-translational energy exchange is 
larger than unity. If the sum of the probabilities is greater than unity, an alternative 
methodology is implemented as shown in Figure 2.2. These selection methodologies are 
easily updated to include electronic energy exchange. 
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Figure 2.2: Selection methodology for inelastic collisions if the sum of the 
probabilities is greater than unity. 

2.2.3 Chemical Reactions 

Once the collision partners have been selected as described in Section 2.2.1, they are then 
tested to see if a chemical reaction will occur. A typical bimolecular reaction may be written 
as: 


A + B C + D 

As long as the reaction takes place in a single step with no species other than the reactants 
present, it has been found that the rate equation for the change of concentration of species A 
may be written: 


--^ = k f (T)n A n B -k r (T)n c n D (2.2.36) 

The rate coefficients are functions of temperature and are independent of the number 
densities and time. A further result that originally had an empirical basis is that the rate 
coefficients are of the form: 
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(2.2.37) 


k f = aT h exp 


kT 


where a and b are constants and E a is called the activation energy of the reaction. 

The Total Collision Energy (TCE) model[55] of Bird was developed in order to obtain a rate 
coefficient based on a chemical reaction cross section which takes the fonn: 
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where Ci and Cj are constants, E co u is the collision energy available to the reaction, and £ is 
the average number of internal degrees of freedom that contribute to the collision energy. 

The resulting rate coefficient becomes: 
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where the constants C/ and C? are: 
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C 2 b 1 + co ab 


The probability of the reaction then occurring at each collision with E c > E a between power- 
law VHS particles is then: 
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The same can be done for termolecular, dissociation/recombination reactions represented as: 

AB + T <-» A+B + T 

where AB is the dissociating molecule and T is the third body particle. The rate of formation 
of either species A or B in the forward reaction can be written: 


dn 


= k f n AB n T 


(2.2.42) 
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while the loss of species A in the reverse reaction is 


dn A , 

— = k r n A n B n T 
dt 


(2.2.43) 


The corresponding reaction probability for the forward dissociation reaction is the same as 
that in Equation (2.2.41) and the probability of a recombination can be written: 
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where the activation energy and, therefore, contributing internal degrees of freedom are taken 
as zero. 

These are the typically implemented set of probabilities for both forward and reverse 
reactions in DSMC codes. They require available data for each forward and reverse reaction. 
However, these data are not always readily available for the reverse reactions. Also, many 
CFD codes make use of statistical mechanics in order to specify the reverse reaction rate. If 
comparisons are to be made to CFD results, it follows that the DSMC code should use 
statistical mechanics to obtain the reverse reaction rate where data are not available. This 
methodology is not used in general where data are available because the reaction 
probabilities are functions of the macroscopic temperature, which is less physically realistic 
than the dependence on the collision energy for the forward reactions. Another use of 
statistical mechanics is to calculate the reverse reaction rate from the forward reaction data, 
and then choose the relevant coefficients from the available reverse reaction data to most 
closely match the rate computed from statistical mechanics. The following modifications 
were implemented in the DAC software for those reverse reactions where data were 
unavailable. 

Chemical equilibrium requires a balance between the forward and reverse reactions such that 
the concentrations of all species remain constant with time. Therefore, for a bimolecular 
reaction, 


- _ k f = n c n D 
^ K n A n B 


(2.2.45) 


where K eq is the equilibrium constant. Similarly, the equilibrium constant for the 
tennolecular dissociation-recombination reaction is written as: 
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(2.2.46) 


Equilibrium statistical mechanics relates the numbers of the various species in a system to the 
partition functions. For the bimolecular reaction, 
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where Q is the partition function for the species that is indicated by the superscript, and E a 
should be interpreted as the difference between the forward and reverse activation energies. 
The number N is related to the volume V of the system by N = nV, which results in the 
bimolecular equilibrium constant written as: 
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For the termolecular reaction, the system volume does not cancel and the corresponding 
equation is: 
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The partition function for each species can be written as the product of the separate 
translational, rotational, vibrational, and electronic partition functions: 


Q = Q tmns Qro'Q vib Qe, (2-2.50) 

First, the translational partition function is: 

Q„„ = v(2mnkT/h 2 f (2.2.51) 

The rotational partition function is: 

<2 2 52) 


where e’ is a symmetry factor equal to 2 for homonuclear molecules and 1 for heteronuclear 
molecules. The harmonic oscillator model leads to a vibrational partition function of: 

e,» = l/(l-«P(-0„ „/T)) (2.2.53) 


Finally, the electronic partition function is summed over the j levels as: 
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where g, are the degeneracies. Finally, it can be shown[2] that, for a bimolecular reverse 
reaction given the data for the forward reaction ( a ’s and />’ s), the reverse reaction probability 
is: 
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Also, for the termolecular reverse reaction: 
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For the case where the reverse reaction has finite activation energy, the probability is then: 
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2.2.4 Electric Field 
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Various researchers have modeled the electric field structure in a hypersonic shock layer in 
the past. The first such calculation was by Bird[26]. In his work, the electrons were moved 
through the computational grid with the ions that they were created with. No explicit electric 
field calculation was made, so the ions and electrons were not accelerated, and the simulation 
time step was that of the heavy particles. Later, Bird modified this method[27] to include a 
calculation of the ambipolar electric field using a form of the equation derived by Langmuir 
and Tonks[56] 


^ ambipolar 
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(2.2.58) 


In this equation, T e and n e are the average electron translational temperature and number 
density, macroscopic quantities derived from the DSMC solution. This equation can be 
derived from the macroscopic conservation of momentum for the electron species using the 
assumption of negligible inertial effects, negligible friction due to collisions, zero magnetic 
field and constant translational temperature. In his work, Bird calculated E am bi P oiar using the 
results from a previously converged calculation. He then re-converged the DSMC 
calculation, applying the prescribed acceleration to the charged particles at each time step 
from E ambipolar ■ The field E am bi P oiar was calculated from the re-converged DMSC solution, and 
the entire process repeated until there was no change in the resulting flow field parameters. 
The movement of the electrons was still tied to that of the ions in this method. 


In the work by Gallis and Harvey [28], the electrons and ions were again moved together, and 
the electric field was calculated in a similar manner to equation (2.2.58). Here, however, the 
assumption of isothermal electrons was not made: 
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and the gradient of the electron pressure was instead used to calculate the ambipolar electric 
field. Gallis and Harvey used a simulation time step corresponding to the cell crossing time 
of the heavy particles, and accelerated the heavy particles according to the field E am bipoiar- 
They then computed the average ion velocity in a given cell, and adjusted the average 
electron velocity in the same cell to match. In this way, the electron energies were affected 
by the electric field in an average sense. 

Carlson, Hassan and Taylor[25, 29] devised a method in which assumptions of zero net 
current and charge neutrality were used to calculate the electric field without the use of any 
macroscopic quantities. They wrote down Newton’s law for a charged particle of species 5 
moving under the influence of an electric field: 

= (2.2.60) 
at 

and thus the expressions for average electron and ion velocity by summing over individual 
particles: 
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where E is the average electric field over a simulation time step, N is the number of 
simulator particles, v 0 is a velocity vector at the start of a simulation iteration, and the 

subscripts e and i refer to electron and ion species, respectively. By enforcing charge 
neutrality, N e = A), they then solved for the ambipolar electric field required to make the 
average ion and electron velocities equal, while enforcing the assumption of zero net current: 
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The electric field given above was used to accelerate the charged particles during the move 
portion of the DSMC algorithm. Additionally, a simple model of the particle acceleration in 
the plasma sheath at the vehicle surface was implemented. The electrons were reflected 
specularly from the vehicle surface to simulate their reflection in the strong negative 
potential gradient that would exist in the plasma sheath. An additional energy increment, e(p, 
was added to the surface heat transfer measurement for each ion impacting the surface. The 
expected potential drop across the sheath, 0, was calculated from one-dimensional 
collisionless sheath theory. 

Boyd[30] developed a model for the electric field of a weakly ionized plasma that is similar 
to the one originally proposed by Bird. The model was first implemented in a simulation of 
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the plasma flow through an arcjet type thruster. Each electron is moved with the velocity of 
an ion throughout the domain, and the charged particles do not receive any velocity 
increment because of their response to an electric field. However, in this model, the average 
ion velocity of the cell in which the electron is located is used to move the electron particles, 
and it is not associated with a specific ion particle. This model maintains charge neutrality in 
an approximate sense, and results in lower computational overhead than the first model 
proposed by Bird. 

In order to calculate the self-induced electric field without any assumption of charge 
neutrality, zero net current, or ambipolar diffusion, the electrostatic Poisson equation must be 
solved. This was carried out for very rarefied flows surrounding spacecraft in low Earth orbit 
by two groups of researchers[57, 58]. While interesting, the preliminary results presented by 
these researchers in the cited references are for flow fields of significantly different structure 
than the shock layers being considered in this work. Gallis et al. [59] presented preliminary 
results in which they solved the electrostatic Poisson equation in conjunction with the DSMC 
method for a hypersonic helium flow field. Unfortunately, the results presented in this work 
were sparse, and the physical processes associated with reacting air were not included in the 
analysis, and a systematic study of the effect of rigorously including the electric field in the 
computation on the flow field parameters was not presented. 

Presently, a hybrid of Bird’s original model and Boyd’s model[30] was implemented. The 
electrons were kept separate but retained their own velocities as was done in Boyd’s model; 
however, they are moved along with their associated ion as was done in Bird’s model. By 
doing this, all particles are moved at the heavy particle time scale, but collisions occur at the 
electron time scale. 

2.3 Recent Direct Simulation Monte Carlo Advances 

Molecular-level chemistry models that predict equilibrium and non-equilibrium reaction rates 
using only kinetic theory and fundamental molecular properties (i.e., no macroscopic 
reaction-rate information) have recently been proposed[ 10, 11, 13, 14], In addition, a new 
vibrational relaxation model has also been introduced that only relies on fundamental 
molecular properties and a single measure of vibrational collision number at a reference 
temperature [10], Preliminary evaluations [14] indicate that the resultant reaction rates are in 
very good agreement with measured Arrhenius rates for near-equilibrium conditions and with 
both measured rates and other theoretical models for far-from-equilibrium conditions. 

2.3.1 Vibrational Relaxation 

A new vibrational relaxation model has recently been introduced by Bird [10], The data for 
the vibrational collision number, Z vi b, is generally presented as a function of temperature in a 
Millikan-White plot. A convenient way of using this information in a DSMC procedure is 
through an expression that employs the values of Z vi b at two temperatures. The temperatures 
may be set to the characteristic temperature of vibration, 0 V7 *, and the characteristic 
temperature of dissociation, ©dm- Then, if the corresponding collision numbers are Z c and 
Zdiss, the variation of the collision number with temperature T is given by: 
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(2.3.1) 
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The problem with the Millikan- White data for most gases is that it leads to collision numbers 
less than unity at the dissociation temperature. This is physically impossible and the most 
effective way to deal with the problem is to set the collision number to unity at the 
dissociation temperature. Also, it is preferable to set the second temperature, T re f, where the 
collision number is Z re f. The collision number then becomes: 
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(2.3.2) 


For a procedure based on this equation, the required data for each vibrational mode is the 
vibrational collision number at the reference temperature. The one-third-power law is 
preserved, but there is a decrease in the rate at higher temperatures. This is essentially a 
more precise alternative to the “Park correction factor ’[60] that is commonly employed in 


CFD. 


Consider a collision in which the sum of the relative translational energy and the pre- 
collision vibrational energy of the molecule under consideration is E co u. The “quantized 
collision temperature” is defined by: 
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(2.3.3) 


where the collision energy has been quantized through the truncation of i max to an integer. 

The relevant parameters for the new vibrational model used are listed in Table A.3. 

2.3.2 Chemical Reactions - The Quantum-Kinetic Chemistry Model 

Bird has proposed[10, 11, 13] a set of molecular-level chemistry models based solely on 
fundamental properties of the two colliding particles: total collision energy, quantized 
vibrational levels, and molecular dissociation energies. These models link chemical-reaction 
cross-sections to the energy-exchange process and the probability of transition between 
vibrational energy states. The Larsen-Borgnakke[42] procedures and the principle of 
microscopic reversibility are then used to derive simple models for recombination reactions 
and for the reverse (exothermic) reactions. These models do not require any macroscopic 
data, and they balance the fluxes into and out of each state, thus satisfying microscopic 
reversibility. As mentioned earlier, preliminary evaluations [14] indicate that the resultant 
reaction rates are in very good agreement with measured Arrhenius rates for near-equilibrium 
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conditions and with both measured rates and other theoretical models for far-from 
equilibrium conditions. 

2. 3. 2 . 1 Dissociation Reactions 

Bird suggests[10] that a complete and consistent description of vibrational energy exchange 
between colliding molecules must include dissociation. This concept is embodied in the 
energy-threshold model of Bird[2]. According to this model, during the energy-exchange 
process of a diatomic molecule, if the energy content of the vibrational mode exceeds the 
energy threshold of the dissociation or exchange reaction, the reaction AB + T — > A + B + T 
occurs. 

In his critique of Bird’s model, Lord[61] pointed out that a “continuum” of states (i.e., an 
infinite number) exists above the dissociation limit. Therefore, Lord argues, if a dissociation 
reaction is energetically possible, it occurs. Thus, assuming a particular collision between 
two simulators, where at least one of them is a diatomic molecule, the serial application of 
the Larsen-Borgnakke energy exchange model would make the collision energy: 

^ coll ~ ^ trails .pair ^ vib , mol (2.3.4) 

available to the vibrational mode of the molecule in question. In this expression of collision 
energy, E tmns ,pair is the translational energy of the pair, and E vi b, mo i is the vibrational energy of 
the molecule considered for dissociation. Assuming a harmonic oscillator model, the 
maximum vibrational level that can be reached is that given in Equation (2.3.3), where the 
vibrational collision number has been defined as unity. If this level is higher than the 
dissociation level (id = ©d / 0,), or 
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(2.3.5) 


a dissociation occurs. 


When the serial application of the quantum Larsen-Borgnakke model is considered, potential 
post-collision states i are chosen uniformly from the states equal to or below i max and are 
selected through an acceptance-rejection routine[2] with probability 


P‘ = 

1 vib 


i*k 0 


vib 


“ coll ) 


(2.3.6) 


where co is the temperature exponent of the coefficient of viscosity. Then, the corresponding 
rate coefficient in an equilibrium VHS gas at temperature T is: 

n) = gC,Wv,r J ' (2.3.7) 

The degeneracy g of the reaction is included to cover the case of a polyatomic molecule in 
which the reaction can occur through more than one mode with the same characteristic 
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vibrational temperature. The collision rate parameter R^ r is the collision rate for collisions 
between gas species AB and T divided by the number density product. It can be written for 
an equilibrium VHS gas as 
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where r re f is the molecular radius at temperature T re f and m r is the reduced mass. For 
dissociations, it is necessary to add the symmetry factor e that is 1 for homonuclear and 2 for 
heteronuclear molecules. The parameter T(i max ) AB J is the fraction of collisions between AB 

and T that have sufficient energy to meet the condition of equation (2.3.5). With i as the pre- 
collision vibrational state of AB, the result for an equilibrium VHS gas is 
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(2.3.9) 


Here, Q(a,x) = F(a,x)/r(a) is a form of the incomplete Gamma function and 

Z.(T)= i/D - exp(-0 v /T )] is the vibrational partition function. As an example, the 

dissociation rate of molecular nitrogen is presented in Figure 2.3 for the Q-K model as 
compared to several rates from the literature, the first of which is the value used in the 
DSMC simulations using the TCE chemistry model. 




Temperature (K) 


Figure 2.3: Comparison of Q-K dissociation rate of the reaction N 2 + N 2 — > N + 
N + N 2 with those from literature. 
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2.3. 2.2 Recombination Reactions 


The condition for recombination in a collision between species A and B is that another atom 
or molecule is within the “ternary collision volume” V co u- The probability of recombination 
is therefore 


Prec = nV coll (2.3.10) 

where n is the number density. The corresponding rate equation is 

K(T) = K;!nV„„ (2.3.11) 

The ratio k f (T)/ k r (T) from Equations (2.3.7) and (2.3.1 1) may be equated to the 

equilibrium constant from statistical mechanics in order to determine the ternary collision 
volume. It has been found that a volume in the fonn 

V co n = aT b V ref (2.3.12) 


allows nearly an exact fit of this ratio to the equilibrium constant. The coefficients used for 
these adjustments can be viewed in Table A. 16. The reference volume F, e /is most 
conveniently set to the volume of the sphere with radius equal to the sum of the reference 
radii of the three particles. 

The values of the constants a and b are determined in the context of the macroscopic 
temperature, but the evaluation of the ternary collision volume in the DSMC implementation 
is based on the quantities associated with the collision. The “collision temperature” based on 
the relative translational energy in a collision with relative speed g is: 
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1 coll 
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(2.3.13) 


As an example, the cumulative recombination rate of two atomic nitrogen atoms is presented 
in Figure 2.4 as compared to values from the literature, the first of which is the rate 
computed by statistical mechanics from the dissociation of molecular nitrogen. 
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Figure 2.4: Comparison of Q-K recombination rate of atomic nitrogen with 
those from literature. 

2.3. 2.3 Exchange Reactions 

Exchange reactions are binary reactions written as A + B C + D in which A and C are the 
molecules that split, while B and D are the atoms. The phenomenological DSMC procedure 
for this reaction is analogous to that for the dissociation model. For both the forward and 
reverse reactions, the reaction probability is equal to that of Equation (2.3.6) with i set to the 
vibrational level corresponding to the activation energy of the reaction, i a 


This is a relative probability that must be normalized through division by the sum of 
probabilities from the ground state to i max , resulting in 


This is sufficient if the level i a is much larger than unity, but some endothermic reactions and 
all exothermic reactions have either zero or small activation energies. The use of discrete 
levels then leads to unacceptably abrupt changes in the rate coefficients and it is impossible 
to obtain a near exact agreement between the ratio of forward to reverse rate coefficients and 
the equilibrium constant. The level i a is the highest level with energy below the activation 
energy E a and a smooth curve is obtained for the rate coefficients if the energy i a kO y in the 
numerator is replaced by E a . The DSMC probability of an exchange reaction in a collision 
with E c > E a is then 


i a = Int 



(2.3.14) 



(2.3.15) 
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(2.3.16) 
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An exact expression can be written down for the fraction of collisions that satisfy the 
condition of Equation (2.3.16) in an equilibrium VHS gas. With j as the pre-collision 
vibrational level, this is 



and the corresponding rate coefficient is 

k(T) = gR cM (dN/N) (2.3.18) 

This is a complex expression not easily programmed, but it is possible to develop an 
alternative expression that is based on physical reasoning rather than formal kinetic theory. 
The selection of possible post-collision vibrational levels in Equation (2.3.6) is from an 
equilibrium distribution so that the rate coefficient is simply the product of the degeneracy g, 
the collision rate parameter, and the fraction of the vibrational level i a in the Boltzmann 
distribution. The term i a k& v in the exponential is replaced by E a in order to maintain 
consistency with Equation (2.3.16), and the normalizing factor is simply the vibrational 
partition function, thus 

k,(T) = «X? ex p(-fy ) (23. W) 

and 

K C T ) = «X;f exp (-^} 4- ( T)‘ (2.3.20) 

The default activation energy for the forward reaction is the heat of reaction and that for the 
reverse reaction is zero. The ratio of the forward to the reverse rate is compared with the 
equilibrium constant and the two are brought into the best possible agreement by the upward 
adjustment of either the forward or reverse activation energy. These adjustments are 
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temperature dependent through a relation similar to that in Equation (2.3.12) for the ternary 
collision factor: 
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(2.3.21) 


For the endothermic reaction, E exc h is the negative of the activation energy, which would then 
result in the second part of the equation adding the prescribed energy to the activation 
energy. For the exothermic reaction, E exc h is equal to the activation energy of the 
endothermic reaction, so the second part of the equation above does not come into effect. 

The coefficients used in these adjustments are presented in Table A. 17. In addition to these 
phenomenological adjustments, there may be experimental or theoretical evidence for 
activation energies above the default values. Should a phenomenological adjustment fall 
short of the desired value, both the forward and reverse activation energies may be adjusted 
to keep the forward to reverse rate ratio in agreement with the equilibrium constant. As an 
example, the endothermic and exothermic reaction rates for the exchange reaction N 2 + O — > 
NO + N are compared to rates from the literature in Figure 2.5 and Figure 2.6, 
respectively. The Q-K results presented in these figures are in good agreement with those 
presented from the literature. 
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Figure 2.5: Comparison of Q-K exchange rate for the endothermic reaction 
N 2 + O — > NO + N with those from literature. 
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Figure 2.6: Comparison of Q-K exchange rate for the exothermic reaction 
NO + N — > N 2 + O with those from literature. 
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Section 3 

Extending the Quantum Kinetic Chemistry Model 


For missions where the shock-layer temperature reaches tens-of-thousands of degrees Kelvin, 
the importance of including electronic energy levels and ionization reactions becomes more 
pronounced. Although the treatment of electronic energy levels and ionization using the 
DSMC method is not a new development, most of these methods are based on measured, 
equilibrium rates, which are always questionable at high temperature, especially when 
applied to non-equilibrium problems. 

In Section 2, the Quantum-Kinetic Chemistry Model of Bird was introduced. One of the 
attractive properties of this model is the fact that it does not rely on any macroscopic 
information, leveraging instead off of kinetic theory and physical data corresponding to the 
molecules undergoing each individual collision. The purpose of this Section is to introduce 
new methods of treating electronic energy level transitions and chemical reactions involving 
charged species following the Q-K methodologies utilizing not only the vibrational energy 
levels, but also the electronic energy levels of the atoms and molecules. 

3.1 Detailed Electronic Energy Level Model 

To model the distribution of electronic energy in the DSMC technique, three procedures must 
be defined. First, when a particle is introduced into a simulation, it is necessary to obtain a 
new electronic energy from a given distribution through statistical sampling. A new 
electronic energy is also required when a particle changes its electronic energy after a 
collision with a surface. Under equilibrium conditions, the distribution has the well- kn own 
Boltzmann form. Second, it is necessary to statistically sample a new electronic energy 
following a collision, which involves electronic energy transfer. Third, a method to 
reproduce the electronic energy transition rate must be defined. Detailed examples of 
applying this model at hypersonic (re)entry conditions in air are presented subsequently in 
Section 4.3. 


3.1.1 Equilibrium Sampling 


Each electronic energy level j has a distinct energy, Ej, and degeneracy, g,. The Boltzmann 
distribution for the electronic energy levels at a given temperature T gives the following 
result for the fraction of particles in level j: 


f el ( 2 ) 


N J= gj^M-Ej/kT) 
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(2.3.22) 


This distribution is used when either creating a new particle or after a surface collision in the 
DSMC simulation at a boundary specified at the temperature T. However, it is not possible 
to sample an electronic energy level j directly from the distribution. Therefore, an 
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acceptance-rejection procedure is employed by selecting values for j from the following 
distribution: 


, () f,ex p(~E,/w) 

gj exp(-Ej/kT) 


(2.3.23) 


where / is the value of / for which Equation (2.3.22) is a maximum. Unfortunately, since the 
degeneracy of each level is a variable specific to each species, there is no way to detennine 
the maximum level a priori, so the maximum level must be searched for at each 
implementation or saved and recalled for a constant boundary temperature. The sampling of 
a new electronic energy level then proceeds as follows: 

1) select at random an electronic energy level evenly distributed between 0 and J max 
where J max is the maximum possible energy level; 

2) detennine the value J for which Equation (2.3.22) is a maximum; 

3) accept the value of j if f e i ’(j) > RAND, a random number; 

4) if the value of j is not accepted, return to step 1 . 

3.1.2 Post-Collision Sampling 


A phenomenological approach is usually adopted in the DSMC method when a collision 
occurs that involves energy transfer. The Larsen-Borgnakke method [42] samples a post- 
collision state from a combined distribution of the translational and electronic collision 
energies of the colliding particles. Based on the approach of Bergemann and Boyd[43] and 
Boyd[51], the Dirac delta function is employed to define the distribution of energies in 
Equation (2.3.22) in the following continuous form: 




g.exp(-E ; /w) 

^max 

X& ex p(-' B .-/* T ) 


^ el 

kT 


7=0 



kT J 


(2.3.24) 


Consideration must now be given to the distribution of translational energy of the colliding 
particles. This distribution is naturally affected by the intennolecular model used since this 
detennines the collision probability. For the present study, the variable hard sphere (VHS) 
collision model of Bird[41] is used. However, it is straightforward to develop the 
formulation for an alternative collision model such as the variable soft sphere of Koura[62]. 
For the VHS model, the distribution of the translational energies is: 
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Using Equations (2.3.24) and (2.3.25), the combined distribution for sampling a post- 
collision electronic energy level j ’ from the total collision energy E co u = E tmns +E e i = 
Etrans ’+E e / ’ is: 
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(2.3.26) 


i=0 


In applying the general Larsen-Borgnakke scheme, it is assumed that local thermodynamic 
equilibrium prevails. Therefore, the temperature T in Equation (2.3.26) is constant. Also, the 
total collision energy is constant, so it is only necessary to perfonn sampling of the post- 
collision state from the following distribution form: 
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(2.3.27) 


Again, an acceptance-rejection procedure is used. The normalized distribution that is 
required is obtained by finding the value of j ’ for which Equation (2.3.27) is maximum, J\ 
which is different for each value of E co u. Therefore, the following distribution is obtained: 
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The procedures are implemented as follows: 

1 . given a pair of particles with total collision energy E co u which undergoes electronic 
energy exchange, determine J’; 

2. detennine the maximum allowable electronic energy level obtainable from E co u, J”; 

3 . take J* to be the smaller of J’ and J’ 

4. as described in the previous section on equilibrium sampling, perfonn an acceptance- 
rejection procedure to sample j ’ from: 


g'if’Ecou) 






(2.3.29) 


3.1.3 Electronic Energy Level Transitions 

A transition from electronic energy level i to level j for species A can be written in the fonn 
of a chemical reaction as A' + B — > A 1 + B . Following the Q-K model, the simplest model is 
to assume that the transition occurs if the electronic energy level of A after a trial Larsen- 
Borgnakke redistribution of the collision energy, E co u = E tra ns + E e i, is j =/. Since the 
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transition is treated as a chemical reaction, there are J max transitions to consider where J max is 
the maximum energy level energetically possible. 

In the DSMC code, the process for determining if an energy level transition occurs, and if so, 
to which transition level, is as follows: 

1 . detennine the value of J max and compute the denominator of Equation (2.3.28); 

2. perform acceptance-rejection according to Equation (2.3.28) until a test level is 
accepted. 

This selection process is analogous to the selection process used in selecting a chemical 
reaction from a list of possible reactions. 

3.2 Chemical Reactions Involving Charged Species 

As outlined in Section 2.3.2, Bird[10, 11, 13], supported by Gallis[14], proposes a set of 
molecular-level chemistry models based solely on fundamental properties of the two 
colliding molecules: total collision energy, quantized vibrational energy levels, and 
molecular dissociation energies. These models link chemical cross-sections to the energy 
exchange process and the probability of transition between vibrational energy states. The 
Larsen-Borgnakke procedures and the principle of microscopic reversibility are then used to 
derive simple models. The models do not require any macroscopic data and they function by 
seeking to balance fluxes into and out of each state, thus satisfying microscopic reversibility. 
Since reactions involving only neutral species have been covered thus far, processes 
involving charged species are presented in the current section. Again, more detailed 
comparisons for these reactions pertaining to hypersonic (re)entry in are be presented in 
Sections 4.4 through 4.7. 

3.2.1 Ionization Reactions 

The procedure for ionization is analogous to that for dissociation. The difference is that 
ionization is linked to the electronic energy level instead of the vibrational energy level as 
with dissociation. For ionization, the collision energy is defined as E co ii=E t rans,pair+Eei e , part icie 
where the first term on the right is the relative translational energy of the colliding pair and 
the second tenn is the electronic energy of the particle being considered for ionization. If the 
collision energy is greater than the ionization energy of the particle, ionization occurs. The 
corresponding rate coefficient in an equilibrium VHS gas at temperature T is: 

k i ; n (T)=R^r(j m S B (2.4.1) 

where the collision rate parameter is identical to that of Equation (2.3.8). The fraction 
of collisions between A and B that have sufficient energy to meet the requirement of 
ionization is then: 
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(2.4.2) 
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where z e i(T) A is the electronic partition function for particle A. 

3.2.2 Associative Ionization Reactions 


The forward (endothermic) associative ionization reaction can be thought of occurring as 
A + B — > AB — > AB + e where AB is an intennediate, excited molecule that may lead to 
the ionized species AB + plus an electron. The first step is similar to the recombination 
procedure described in Ref. [10]. If the AB molecule is in the ground vibrational state after 
a trial Larsen-Borgnakke redistribution of the collision energy, the intennediate molecule is 
“formed.” In the second step, the dissociation energy is first added to the original collision 
energy. The new collision energy is then distributed among the internal energy modes of the 
AB molecule and then this molecule is tested for ionization as described above. If the 
remaining collision energy after redistribution plus the electronic energy of the molecule is 
greater than the ionization energy of the molecule, the associative ionization reaction occurs. 
The corresponding rate coefficient can then be written as: 

K=R‘^%,lzATf (2.4.3) 


where the collision rate parameter is again identical to that of Equation (2.3.8). The 
vibrational partition function is included to account for the probability that the AB* molecule 
will be in the ground state. The fraction of collisions with sufficient energy for the molecule 
AB to ionize are then calculated by summing over all possible initial internal energy states of 
atoms A and B (electronic energy states) and all possible vibrational, rotational, and 
electronic energy states of the intermediate AB molecule as follows: 
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(2.4.4) 


The sums over the atomic electronic energy levels is over all possible levels. For the 
molecular internal energies, the maximum level is that which is the highest energetically 
possible. To start, an initial translational energy is assumed to be the average relative 
translational energy, or ( 5 A ~ co)kT . This is an approximation and may result in a 

discrepancy between the above analytical expression and the actual sampled rate. It is 
possible to integrate over all possible relative translational energies, but it makes an already 
long calculation (summing over all possible internal states) much longer. As will be shown 
subsequently, this formulation results in an acceptable comparison with the sampled values 
in Section 4.4. The first two expressions after the summations are simply the probabilities of 
atoms A and B being in electronic levels i e 7 and i e i e B , respectively. The next three tenns are 
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analogous to Equation (2.3.29) in that they are the probabilities of the respective internal 
mode of energy being at the specified quantum level at the specified collision energy. The 
collision energies are defined as: 
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where each successive collision energy subtracts the previous internal energy mode. Then 
the last term, as in the dissociation and ionization rate equations, is the fraction of collisions 
for which the translational energy is greater than the difference between the ionization energy 
of AB and the electronic energy of AB . 

The reverse (exothennic) associative ionization reaction is similar to the forward reaction. If 
the AB * molecule is in the ground electronic state after a trial Larsen-Borgnakke 
redistribution of the collision energy, an intennediate molecule is “formed.” In the second 
step, the ionization energy is first added to the original collision energy. The new collision 
energy is again distributed among the internal energy modes of the AB molecule, which is 
then tested for dissociation as described above. If the remaining collision energy after 
redistribution plus the vibrational energy of the molecule is greater than the dissociation 
energy of the molecule, the reverse associative ionization reaction occurs. The 
corresponding rate coefficient can then be written as: 
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where the collision rate parameter is again identical to that of Equation (2.3.8). The 
last tenn is included to account for the probability that the AB* molecule will be in the 
ground electronic state. The fraction of collisions with sufficient energy for the molecule 
AB to ionize are then calculated by summing over all possible initial internal energy states of 
AB + (the electron has none) and all possible vibrational, rotational, and electronic energy 
states of the intermediate AB molecule as follows: 
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The sums over the internal energy levels of AB' are over all possible levels. For the internal 
energies of AB , the maximum level is that which is the highest energetically possible. In 
order to do this, an initial translational energy is again assumed to be the average relative 


42 



translational energy, or ( 5 A~ co)kT , along with the aforementioned cautions. The first three 

expressions after the summations are simply the probabilities of the AB 1 molecule being in 
the respective states. The next three terms are again analogous to Equation (2.3.29) in that 
they are the probabilities of the respective internal mode of energy being at the specified 
quantum level at the specified collision energy. The collision energies are defined as: 
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where each successive collision energy subtracts the previous internal energy mode. Then 
the last term, as in the dissociation and ionization rate equations, is the fraction of collisions 
for which the translational energy is greater than the difference between the dissociation 
energy of AB and the vibrational energy of AB 

3.2.3 Charge Exchange Reactions 


The charge exchange reaction is analogous to the exchange reaction outlined in Section 
2. 3. 2. 3. The binary exchange reaction can be written as A + B + A + + B . This reaction 
takes place when the electronic energy level of the colliding neutral particle, A, after a trial 
Larsen-Borgnakke redistribution of the sum of the relative translational energy and electronic 
energy of the neutral particle is equal to the electronic energy level with energy just below 
the activation energy of the reaction. The probability of the reaction occurring is then: 


P = - 

chex i m 


o ( e - E a Y 2 

baX^coll ^el) 

..iax 3 

Is, 


1=0 


The corresponding rate equation for the forward reaction is then: 

g.exp (-EF'IkT) 


k chex (T) = R 


coll j A 


&x p(- E 'f/ kT ) 


i = 0 


and for the reverse reaction is: 
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As with the exchange reactions in Section 2.3.2 .3, the ratio of the forward and reverse 
reactions can be compared with the equilibrium constant and adjustments can be made to the 
phenomenological model to bring the model in closer agreement with statistical mechanics. 

3.2.4 Other Exchange Reactions 

In other exchange reactions involving charged species, the predominant mechanism is the 
breaking of a bond of one species (either neutral or charged) in order to exchange one atom 
to bond with the original atom, just as with an exchange reaction involving only neutral 
particles. Therefore, the logic for performing this reaction involving charged species is 
similar to that involving only neutral species. 

3.3 Selecting a Reaction 

Next, the Quantum-Kinetic chemistry model requires the selection of which reaction occurs. 
For many collision pairs, there may be more than one possible reaction. In the TCE 
chemistry model, the reaction probability is computed for each of the possible reactions and 
summed. This sum is then compared to a random number between zero and one. If the 
random number is less than the sum of the probabilities, a reaction is to take place. The 
selection of the reaction then proceeds by normalizing the reaction probabilities, and another 
random number is chosen to select which reaction is allowed. In the Q-K method, there is no 
clear analogy, and previous implementations have utilized a hierarchical approach where the 
reactions are considered by type first. For the neutral reactions, the dissociation reaction was 
always considered first and allowed to proceed even if an exchange reaction was a 
possibility. However, it was found that this approach tended to be biased in favor of 
dissociation reactions, especially at high temperatures where most of the collisions have 
larger relative translational energies and the vibrational level distribution is shifted to the 
higher levels. This method becomes even more undesirable when charged species reactions 
are included in the reaction set. Specifically, what happens when both dissociation and 
ionization reactions occur simultaneously? This event doesn’t actually occur for the current 
reaction set under consideration (only electron impact ionization of atomic species is 
considered herein), but this is obviously an issue that needs to be resolved. 

An alternate approach to the hierarchical method was therefore implemented. This approach 
is, like the TCE model, based on assigning reaction probabilities to all of the possible 
reactions, summing them, then comparing the sum to a random number. If the random 
number is less than the sum of the probabilities, then the probabilities are normalized and 
then another random number is generated to select the reaction. The reaction probabilities 
for exchange and charge exchange reactions have already been given in Equations (2.3.16) 
and (2.4.10), respectively. For the remaining reactions that occur when the collision energy 
is such that the reaction occurs (dissociation, ionization, and associative ionization), the 
probability is set equal to one if it is detennined that the reaction is allowed to occur. These 
probabilities are then summed and compared to a random number between zero and one 
(obviously, if a dissociation, ionization or associative ionization reaction are possible, a 
reaction will occur). If a reaction is to occur, P r i is defined as the reaction probability for the 
i th reaction, and the cumulative distribution for all reactions up to and including the j' h 
reaction is: 
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(2.5.1) 



where NR is the total number of reactions possible with the current reactants. A random 
number is then chosen to determine which reaction proceeds. 

As shown subsequently in Section 4, this selection methodology allows the sampled reaction 
rate to accurately reproduce the analytical rates shown in this Section. There is one item that 
requires mentioning at this point that the chemistry set utilized in this thesis does not include. 
That is the possibility of a dissociation and ionization reaction occurring for the same set of 
particles. For example, when molecular nitrogen collides with an electron, the molecule can 
either dissociate or ionize. Under the rules for these processes listed above, when the 
respective criterion is met for either reaction, the probability of that reaction occurring is 
unity. What happens when both have a probability of unity? Using the above reaction 
selection methodology, it is shown in Figure 3.1 that the sampled rates still accurately 
reproduce the analytical rates. This occurs because the energy required for dissociation is 
less than that for ionization, so will always occur more frequently as shown in the figure. 

The highest reaction rate physically possible would be equal to the collision rate of the 
particles. As the limit where all of the collisions have sufficient energy for both reactions is 
approached, a departure from the analytical rates occurs since both are equal, resulting in 
sampling rates one half the analytic rates. 




Figure 3.1: Sampled and analytic reaction rates when dissociation competes 
with ionization. 
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Section 4 

Application of New Models to Upper-Atmosphere Hypersonic Flows 


In this Section, the new electronic energy level transition and chemistry models are validated 
by examining their ability to reproduce rates of chemical reactions in the literature [5 -7, 9, 15, 
29, 60, 63-78]. The DSMC models are applied in a two-dimensional code where the results 
are accumulated over all cells to approximate a zero-dimensional simulation. The test gas is 
comprised of eleven species: O, N, O 2 , N 2 , NO, 0 + , N + , 0 2 + , N 2 + , NO+, and e'. The 
simulations were run at various temperatures between 10,000 K and 60,000 K and at a 
continuum number density of 1.0e23 /nr in an adiabatic box 0.002 m on a side with 2 million 
particles. During the simulations, relaxation of the particles (including rotational, vibrational 
and electronic) was allowed to proceed as usual, but when a reaction was determined to take 
place, the number of reactions was advanced by one but the simulators were left unchanged, 
to maintain the thermodynamic energy balance with energy neither added nor removed from 
the gas. 


4.1 Measurement of Reaction/Transition Rates 

For a bi-molecular reaction A + 5<->C + h, the forward reaction rate is calculated as: 


dn , 


k f = - 


dt 


(3.1.1) 


where the numerator is the change in number density of species A due to the reaction in the 
sampled time and the denominator is the product of the number densities of species A and B. 
The corresponding reverse reaction rate is: 


dn. 


k, = ■ 


dt 


n c n D 


(3.1.2) 


For a termolecular reaction AB + T <-> A + B + T , the forward reaction rate is calculated as: 

dn,„Z 


k f = 


dt 


(3.1.3) 


as the bi-molecular reaction rate is, but the reverse reaction rate is now: 

dn. / 


k, = ■ 


dt 


2 n A n B n T 


(3.1.4) 


The capability of these techniques to reproduce measured equilibrium reaction rates will now 
be investigated. However, the uncertainty of the currently available measured rates for 
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atmospheric hypersonic flows can be large. Most of the experimental data on which reaction 
rate sets are based were obtained for low temperatures (2,000 - 7,000 K). Even though 
uncertainty may be large, Park[60] suggests that the experimentally obtained reaction rates 
can be extrapolated up to 30,000 K. 

The uncertainty of these measured rates is not always reported and in some cases is 
unknown. The problem becomes even more severe in the extrapolated portion of the 
temperature range (T> 10,000 K). Park[60] reports that some rates measured for the same 
reaction differ by more than one order of magnitude, which will be apparent in the figures 
that follow. Since this work is focused on the extension of the Q-K model to electronic 
energy levels and reactions involving charged particles, only these transition and reaction 
rates will be discussed in this Section. However, comparisons of the analytical Q-K rate 
expressions, the sampled Q-K rates, and rates found throughout the literature for every 
reaction included in this study are presented in Appendix C. 

4.2 Definition of Translational, Rotational, Vibrational, and Overall 
Temperatures 

There are various ways to define temperatures with regards to the DSMC methodology. In 
this work, we will refer to statistical mechanics for the definition of temperatures by defining 
the energy as: 


e = %£KT (3.2.1) 

where £ is the specific energy of the quantity of interest, £ is the number of degrees-of- 
freedom, R is the universal gas constant, and T is the temperature. By finding the energy, 
recognizing that k = mR , and rearranging Equation (3.2.1), one has the temperature of 
interest. For the translational temperature, we define e t rans as: 

= ( 3 - 2 - 2 ) 

where c' 2 is the average square of the thennal velocity defined as: 

c' 2 = c 2 - 2cc 0 + c] (3.2.3) 

where c is the mean molecular velocity and c 0 is the macroscopic stream or mean velocity. 
The translational temperature can then be defined as: 

2e 

T = ,rans n 7 41 

1 trans 

since there are 3 translational degrees of freedom. 

For the internal degrees of freedom of a specific species n, we define the average internal 
energy of a single particle as: 
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(3.2.5) 


77’ n _ n 

^int = m n £ m 


and the temperature is then: 
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int b-K n 
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(3.2.6) 


For the overall internal temperature averaged over all of the N species, we then define: 

N 

\ ’ fn rjin 

/ j b inf* int 


T 


Zfi 


(3.2.7) 


For rotational energy, the average energy of each species is found by summing the rotational 
energies of the particles in the computational cell and dividing by the number of particles. 

For the rotational degrees of freedom, it is assumed that the mode is fully excited and is set to 
2 for diatomic molecules. 

For vibrational energy, the average energy of each species is again found by summing the 
vibrational energies of the particles and dividing by the number of particles. To then find the 
number of vibrational degrees of freedom for the species, using the definition of the average 

vibrational energy level i n = E" ib /kQ vib and some manipulation, one gets: 

C. = 2fin(l + 1 II) (3.2.8) 


The electronic temperature is treated in a similar manner, but some modifications must be 
made to account for the discretized spacing of the electronic energy levels. To begin, the 
species electronic temperature is defined as: 
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ele ^ ele 


1 ele 


k In I 


«o£'i 

i n n 

V n xS 0 


(3.2.9) 


where the 0 and 1 denote the ground and first excited electronic energy levels, g is the 
degeneracy of the level, and n can be taken as the number density or total number of samples 
(result is the same since the ratio is being taken). This definition comes directly from the 
ratio of the Boltzmann distribution of the ground and first excited electronic energy levels. 
Then, the effective number of degrees of freedom is defined as: 
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as we have defined the effective degrees of freedom of the other internal modes. 

4.3 Electronic Energy Level Transitions and Comparison to 
Equilibrium 

The following section applies the methodologies developed in Section 3.1 to the 
aforementioned equilibrium heat bath. There are excitation and quenching processes 
between each of the atomic/molecular energy levels, resulting in a large number of rates 
being required for previous electronic energy level transition models, both for CFD and 
DSMC. The new Q-K based transition model negates the necessity for each of these 
transition rates (most of which use approximate analytic formulas). In the following sub- 
sections, only a portion of the transitions presented in Appendix C are discussed for 
exemplary purposes only. 

4.3.1 Electron-Impact Excitation 

The process of electron-impact excitation is significant for both atoms and molecules in 
hypersonic shock-layer applications. To begin, the electron- impact transitions of atomic 
nitrogen are examined. The sampled transition rates are compared to values reported by 
Tayal[79], Park[60], and Frost[80]. The transitions from the first to the fifth level and from 
the second to the fifth level are presented in Figure 4.1 and Figure 4.2, respectively. It is 
important to note that the 1-5 transition is considered a “forbidden” transition, while the 2-5 
transition is an “allowed” transition. Even so, the phenomenological Q-K model compares 
very well to the reported rates from the literature. 




Temperature (K) 


Figure 4.1 : Comparison of sampled and reported electronic energy level 
transition rates from the first to fifth energy levels of atomic 
nitrogen. 
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Figure 4.2: Comparison of sampled and reported electronic energy level 
transition rates from the second to the fifth energy levels of 
atomic nitrogen. 

For the electron-impact transitions of atomic oxygen, the sampled transition rates are 
compared to values reported by Zatsarinny[81], Park[60], and Bhatia[82]. Once again, even 
though one transition is allowed (level 1-5, Figure 4.3) and one is forbidden (level 2-5, 
Figure 4.4), the comparisons are quite favorable. 



Temperature (K) 


Figure 4.3: Comparison of sampled and reported electronic energy level 
transition rates from the first to fifth energy levels of atomic 
oxygen. 
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Figure 4.4: Comparison of sampled and reported electronic energy level 
transition rates from the second to the fifth energy levels of 
atomic oxygen. 

The electron-impact transition rates of the molecular species are compared to rates presented 
by Chemi[83], Telulet[84], and Gorelov[85]. Example transition rates for molecular 
nitrogen, nitric oxide, and the molecular nitrogen ion are shown in Figure 4.5, Figure 4.6, 
and Figure 4.7, respectively. Again, good agreement is shown between the sampled rates 
and the reported rates. 




Temperature (K) 


Figure 4.5: Comparison of sampled and reported electronic energy level 
transition rates from the second to the third energy levels of 
molecular nitrogen. 
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Figure 4.6: Comparison of sampled and reported electronic energy level 

transition rates from the first to the fourth energy levels of nitric 
oxide. 



Figure 4.7: Comparison of sampled and reported electronic energy level 
transition rates from the first to the third energy levels of the 
molecular nitrogen ion. 

4.3.2 Heavy Particle Impact Excitation 

Transitions due to heavy particle collisions are similar to electron-impact transitions, but they 
are significantly less important for most applications, as discussed by Johnston[86], where he 
shows that the process is negligible for atomic and molecular species in lunar-return shock 
layers. The comparisons shown in this section are, for the most part, as favorable as those in 
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the previous section for electron-impact transition rates. Comparisons are made with values 
reported in References [85, 87-93]. For example, the transitions of the molecular nitrogen 
ion when colliding with molecular nitrogen, as shown in Figure 4.8, is generally within an 
order of magnitude of the reported values (note the 3-4 order of magnitude in the spread of 
data for the level 1-2 transition). There is also good agreement for collisions of molecular 
nitrogen and nitric oxide for most transitions, as exemplified in Figure 4.9 and Figure 
4.10, respectively. 




Temperature (K) 


Figure 4.8: Comparison of sampled and reported electronic energy level 

transition rates for the molecular nitrogen ion when colliding with 
molecular nitrogen. 
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Figure 4.9: Comparison of sampled and reported electronic energy level 

transition rates for molecular nitrogen when colliding with atomic 
nitrogen. 



Figure 4.10: Comparison of sampled and reported electronic energy level 
transition rates for nitric oxide when colliding with molecular 
nitrogen. 

However, there seem to be certain combinations of heavy-impact partners where the current 
phenomenological model can’t accurately reproduce the measured/calculated transition rates. 
These combinations generally include transitions where the collision partner is molecular 
nitrogen, as shown in Figure 4.11, and the transition is from the lower levels (i.e., 1-2 or 
even some times 2-3). 
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Figure 4.11: Comparison of sampled and reported electronic energy level 

transition rates for the nitric oxide when colliding with molecular 
nitrogen. 

The current model also appears to have some difficulty in reproducing electronic energy 
level transitions in collisions of molecular oxygen, as shown in Figure 4.12. However, 
even though these comparisons have been less than ideal, they are not completely unreliable 
considering the four to five order of magnitude spread in the rates previously reported from 
the literature. 



Temperature (K) 


Figure 4.12: Comparison of sampled and reported electronic energy level 

transition rates for molecular oxygen when colliding with atomic 
oxygen. 
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4.3.3 Equilibrium Conditions 


In order for the electronic energy level procedures to be acceptable, they need to 
reproduce equilibrium conditions after the simulation has been allowed to run for 
a sufficient amount of time to ensure that detailed balance has been enforced. The 
first requirement is for the sampled translational and electronic temperatures to 
remain very near the input equilibrium temperature. The sampled translational, 
electron, and electronic temperatures for all of the species are presented in 

Table B.1 at 10,000, 20,000, and 30,000 K. All of the temperatures depart from the set 
equilibrium value by less than 0.5%. 

The resultant population distribution for the first seven energy levels of all atomic and 
molecular species (where applicable) for a range of temperatures is compared to the 
Boltzmann distribution in Table B.2 through Table B.31 . Again, the simulation represents 
the equilibrium state very well, with sampled departures from the Boltzmann distribution less 
than 1% for most levels where the number fraction in the state is above 1.0x1 0' 3 . The larger 
departure from the equilibrium distribution is expected for levels with a very small 
population due to statistical noise. 


4.4 Associative Ionization Reactions 

Because of their low threshold energies, associative ionization reactions play an important 
role in determining the level of ionization and structure of the flow field at hypersonic 
reentry conditions. There exists very little experimental or computational data in the 
literature regarding the associative ionization reactions in air. The associative ionization 
reaction N + O NO 1 + e has the lowest energy threshold, however the associative 
ionization reactions of N + N <-» Af + e~ and O + OftOl + e" are also very important 
mechanisms in the current application. The sampled DSMC reaction rates of endothermic 
and exothermic reactions of the N + O reaction are compared to measured rates from the 
literature in Figure 4.13 and Figure 4.14, respectively. 
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Figure 4.13: Endothermic associative reaction rates for A + O — > NO + + e 
compared to rates from the literature. 
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Figure 4.14: Exothermic associative reaction rates for NO + + e N + 0 
compared to rates from the literature. 

In these figures, the black lines are the analytical solutions from Equations (2.4.3) and (2.4.6) 
, the black symbols are the sampled rates from the DSMC code, and the first of the values 
from literature is the curve used by LAURA. As discussed earlier, there is a wide spread in 
the reaction rates reported in the literature. Given the assumed measurement uncertainty, the 
Q-K reaction rates are in good agreement with the rates from the literature. 
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4.5 Ionization Reactions 


Another important reaction mechanism for ionized, hypersonic reentry problems is the 
ionization of an atom or molecule by direct impact. There is very limited data in the 
literature for heavy particle impact ionization and these reactions are not generally included 
in reentry computations, so they will not be discussed in this work. There is, however, an 
appreciable amount of data related to electron impact ionization reactions, especially for 
atomic species. The electron impact ionization reactions generally considered in reentry 
computations are those for atomic nitrogen and oxygen. Reaction rate curves for the electron 
impact ionization reaction for atomic oxygen are presented in Figure 4.15. 
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Figure 4.15: Endothermic ionization reaction rates for 0 + e -> 0 + e + e~ 
compared to rates from the literature. 

4.6 Charge Exchange Reactions 

In charge exchange reactions, a singly charged ion (such as N + or O ) collides with a neutral 
atom or molecule and captures one of its electrons, thereby becoming a neutral atom. Only a 
small amount of energy is transferred to the electron donor, so the newly created neutral 
particle retains most of the original energy of the ion. These can therefore be an important 
mechanism for the creation of high-energy neutral particles. The measured reaction rates 
from the endothermic and exothermic charge exchange reactions corresponding to 
N 2 + N + — > N 2 + N are shown in Figure 4.16 and Figure 4.17, respectively, compared to 
the values from the literature as shown. 
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Figure 4.16: Endothermic charge exchange reaction rates for 
N + + N 2 N 2 + N to rates from the literature. 
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Figure 4.17: Exothermic charge exchange reaction rates for N 2 + N ^ N + + N 2 
compared to rates from the literature. 

4.7 Exchange Reactions Involving Charged Species 

The remainder of relevant chemical reactions involves various exchange reactions involving 
charged species. As an example, the endothermic and exothermic reaction rates for the 
reaction NO + + N <-» N 2 + O are shown in Figure 4.18 and Figure 4.19, respectively, 
compared to values from the literature as shown. 
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Figure 4.18: Endothermic reaction rate for the exchange reaction involving 
charged species NO + + N <-» + 0 with comparisons to rates 

from the literature. 
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Figure 4.19: Exothermic reaction rate for the exchange reaction involving 
charged species N^ + O <-> NO + +N with comparisons to rates 
from the literature. 


60 


Section 5 

Modeling of the FIRE II Reentry Flow with Ionization and 
Electronic Energy Levels with DSMC 

In this Section, the newly developed DSMC models presented in previous Sections are 
exercised on a demanding near-continuum reenty flow field with ionization. These results 
are compared with results from traditional DSMC models, a solution of the FIRE II flow 
field obtained using current best practices in computational fluid dynamics (CFD), and 
values inferred from flight. Both flow field and surface parameters are examined. To this 
end, three DSMC cases are summarized in Table 5.1 . 

Table 5.1: List of DSMC FIRE II cases. 


Case 

Electronic Energy Model 

Chemistry Model 

Case 1 

None 

TCE 

Case 2 

Present 

TCE 

Case 3 

Present 

Q-K 


5.1 The Project FIRE Flight Experiment 

The Project FIRE (Flight Investigation Reentry Environment) was an Apollo-era flight 
experiment to measure the radiative and convective heating on spacecraft material during 
atmospheric reentry at lunar return speeds (1 1 .35 km/s)[2 1]. Project FIRE consisted of two 
sub-orbital flight experiments launched April 14, 1964 and May 22, 1965 from Cape 
Kennedy, Florida. For the current study, computations are performed for the reentry 
conditions of the second flight at 1634 seconds, which corresponds to an altitude of 76.42 
km, and conditions for the simulation are provided in Table 5.2. The geometry and 
dimensions of the second vehicle in the experimental series, FIRE II, are presented in Figure 
5.1 . The forebody of the FIRE II reentry vehicle consisted of three phenolic-asbestos heat 
shields sandwiched between beryllium calorimeters. The first two heat shield and 
calorimeter packages were designed to be ejected after the onset of melting, yielding heating 
data free of the effects of ablation and three separate data gathering periods. Calorimeter 
plugs and radiometers were located at various positions on the heat shields. The flight 
condition examined in this study occurred during the first data collection period. 


Table 5.2: Free stream and surface parameters for FIRE II entry condition. 


Time (s) 

Altitude (km) 

P~ (kg/m 3 ) 


\/„ (m/s) 

M K) 

Kn D 

T wa „ (K) 

Xn2 

Xo2 

1634 

76.42 

3.72E-05 

7.73e20 

11360 

195 

2.59E-03 

615 

0.76 

0.24 


61 

























R 0.9347 m 



Figure 5.1 : Schematic of FIRE II test vehicle. 

5.2 Details of the Continuum Flow Model and Non-Equilibrium 
Radiation Model 

Two different numerical approaches are applied: a CFD solution of the Navier-Stokes 
equations using the Langley Aerothermodynamic Upwind Relaxation Algorithm code and 
particle-based DSMC computations using the DAC code, which is the focus of this study. 

All computations presented herein were steady state and axisymmetric with the x-direction 
along the axis. In each case, an 1 1 -species, 21 -reaction, thennochemical non-equilibrium 
approach is adopted. The surface catalytic boundary condition used in the CFD solution was 
“fully-catalytic,” which means that all atomic species recombine to their molecular 
counterpart (i.e., N + N — > N 2 ) and also charged heavy species recombine with electrons to 

form their charge neutral counterpart (i.e., N + + e~ — » N ). 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is a high fidelity, 
structured grid computer code, specialized for hypersonic re-entry physics, utilizing state-of- 
the-art algorithms for computational fluid dynamic simulations[94, 95]. Key elements of 
LAURA include Roe’s averaging[96] and Yee’s Symmetric Total Variation Diminishing 
(STVD)[97] fonnulation of second-order, inviscid flux. Yee’s STVD fonnulation has been 
found to be exceptionally robust and Courant number independent using first point-implicit 
and then line-implicit relaxation for hypersonic flow simulations. 

A two-temperature thennochemical non-equilibrium model[95] is applied. The chemical 
reaction rates are compiled from previous studies (listed as the references for the first 
reaction rate parameters for each reaction in Table A. 14). The thennophysical properties 
are taken from the work of Mcbride et a/.[98] and Gupta et a/.[76]. Multicomponent 
diffusion is approximated using Sutton and Gnoffo’s[99] approximate-corrected approach. 
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The shock-layer radiation is modeled with the HARA (High-temperature 
Aerothennodynamic RAdiation) code. The details of this code for treating air species are 
presented by Johnson et al. [19, 20]. Briefly, it is based on a set of atomic levels and lines 
obtained from the National Institute of Standards and Technology (NIST) online 
database[100] and the Opacity Project[101], as well as atomic bound-free cross sections from 
the TOPbase[102]. The negative nitrogen and oxygen ions are treated using cross sections 
suggested by Soon and Kunc[103] and Chauveau et al. [104], respectively. The molecular 
band systems are treated using a smeared-rotational band (SRB) model[105], which was 
shown by Johnston et al.[ 19] to be sufficient for treating VUV absorbing and optically-thin 
emitting band systems in air. The molecular data for modeling these band systems are 
obtained from Laux[106], except for the VUV N 2 systems, which are obtained from various 
other sources [107- 109]. The non-Boltzmann modeling of the atomic and molecular 
electronic states is based on a set of electron-impact excitation rates compiled from the 
literature and presented in detail by Johnston et al.[ 20], Following the work of Park[l 10], 
the quasi-steady state assumption is made when solving the Master Equation. The tangent- 
slab approximation is applied to calculate radiative flux and the divergence of the radiative 
flux. 


5.3 DSMC Simulation Parameters 

In this section, the parameters used in the simulation of the FIRE II flight experiment are 
discussed. 

5.3.1 Collision Modeling 

The baseline collision model implemented in DAC, for both neutral and charged species 
interactions, is that originally published by Bird[41], in which all particle interactions are 
treated using the VHS collision model. The reference diameters of the heavy ions 0? + , V? + , 
0 + , N~, and NO are set equal to the diameters of their associated neutral particles, and the 
reference diameter of the electrons is set to le-10 m. This value of the electron reference 
diameter is much larger than the classical value of 6e-15 m, and is used to produce electron 
collision cross sections of the required magnitude within the mathematical framework of the 
VHS collision model. A complete list of the VHS reference diameters for all species and the 
parameters Tref and (O is given in Table A. 1 . As shown previously, the present choice of 
electron reference diameter does not introduce significant error in the DSMC results. 

5.3.2 Internal Energy Modeling 

For the current thesis, the rotational collision number is set to unity. The reason for this is 
because the solutions are being compared to a computational fluid dynamics (CFD) solution 
where the translational and rotational temperatures are assumed equal. Therefore, for DSMC 
to best match CFD, the rotational energy is sampled from the equilibrium distribution at each 
collision (Z rot = 1.0). The vibrational energy relaxation model used is that based on the work 
of Millikan and White and the parameters used for each of the molecular species are 
presented in Table A. 2. 
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One goal of this work is to examine the effects of adding an electronic energy transition 
algorithm to DSMC solutions of hypersonic (re)entry flows. To that end, the electronic 
energy levels of all relevant simulated species (O 2 , N 2 , O, N, NO, 0 2 + , N? + , 0 + , N + , and NO + ) 
are included in Cases 2 and 3. The energy levels included can be viewed in Table A.4 
through Table A.13. The energy levels for O and N were taken from Johnston[86], while 
the remaining species were taken from Hash[l 1 1]. A more comprehensive set of energy 
levels could be compiled for 0 + and N + , but these are the energy levels used by the CFD code 
to be compared with, as well as the nonequilibrium radiation code. 

5.3.3 Chemistry Modeling 

The baseline chemistry model is the TCE model using the reaction rates used in LAURA. 
These rates are listed in Table A. 14 as the first rate listed for each reaction. As discussed in 
Section 2.2.3, the exothermic (reverse) reaction rates are then determined by considering the 
statistical mechanics equilibrium constant, and then the coefficients are chosen from those in 
the literature to most closely match the rate as computed by statistical mechanics. 

The other chemistry model, which is one of the main contributions, is the Q-K chemistry 
model and the extensions to charged species. The characteristic temperatures and reference 
parameters for vibrational relaxation can be found in Table A.1 and Table A.3, 
respectively. Additionally, the parameters required to account for the statistical mechanics 
equilibrium constant and adjust the activation energies for recombination and 
exchange/charge exchange can be found in Table A.16 and Table A. 17, respectively. 
Comparisons of the Q-K reaction set to various rates found in the literature can be viewed in 
Appendix C. 

5.3.4 Electric Field Modeling 

The electric field model implemented in the current study, which was also added to DAC, is 
a hybrid of Bird’s model[26] and that of Boyd[30]. In this model, the electron particles are 
moved along with their associated ion during the movement portion of the DSMC algorithm. 
However, the electrons retain their individual velocity components, which are used in the 
collision portion of the DSMC algorithm. Neither the electrons nor ions actually experience 
an electric field, in that their velocity components are not adjusted to account for an electric 
field as one is not computed. 

5.3.5 Gas-Surface Interactions 

The surface catalycity must be treated differently with the present DSMC simulations 
because the CFD solutions use the fully-catalytic surface boundary condition to account for 
the higher surface temperatures in the continuum regime. However, there is evidence[25] 
that, for surface temperatures such as those at the higher altitudes for the FIRE II experiment, 
the catalycity of neutral atomic species recombining is much closer to zero. Therefore, the 
interaction of particles with the vehicle surface is managed by assuming both a fully catalytic 
surface as was used with the CFD solution for comparison and also a partially catalytic 
surface that only allows charged species to become neutralized. The fully catalytic 
simulations will be referred to simply by their case number (i.e., Case 1) and the partially 
catalytic simulations will have an extra notation indicating that it is partially catalytic (i.e., 
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Case l pc ). The remaining collisions (neutral molecules), along with the recombined particles, 
are then allowed to reflect diffusely from the surface, after thermally accommodating to the 
specified wall temperature. 

5.4 FIRE II Simulations 

In this section, a progression of DSMC solutions is discussed and compared with CFD 
beginning with a traditional DSMC solution without electronic energy levels (Case 1), 
followed by the addition of electronic energy levels (Case 2), concluding with the Q-K 
chemistry model (Case 3). At the end of the section, a discussion of the surface heating, both 
convective and radiative, is included where the DSMC rates are compared to both the CFD 
solution and values inferred from flight. 

5.4.1 Baseline Solution 

The baseline solution, referred to as Case 1 in Table 5.1, is the FIRE II simulation at 1634 
seconds using the TCE chemistry model without the addition of electronic energy levels. 
Various temperatures for this case are compared to that from the CFD solution in Figure 
5.2. Clearly, the exclusion of electronic energy levels from the DSMC solution changes the 
shock-layer profile. The shock standoff distance of the DSMC solution is nearly 40% greater 
than that of the CFD solution, while the shock-layer temperature is also greater than that 
predicted by CFD. It can be seen that, inside the shock layer of the DSMC solution, the 
temperatures are not in equilibrium. This has been reported in Refs. [25] and [4] as an effect 
of the addition of ionization physics at this altitude. There is some scatter present in the data 
for both the rotational and vibrational temperatures, which may be caused by the low number 
density of molecular species in the flow field. The spike in vibrational temperature outside 
of the shock layer is the result of trace species of vibrationally excited molecules that have 
escaped from the shock layer. This is a phenomenon typical of flows such as this in DSMC 
simulations. Since the physics are clearly not the same, comparisons of species 
concentrations are not provided for this case. 
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Figure 5.2: Translational temperature comparison between CFD and Case 1 
(no electronic energy levels). 

5.4.2 Addition of Electronic Energy Levels 


Electronic energy levels were included in DSMC Case 2 and are more appropriate for 
comparison to the CFD solution. The various temperature profiles from CFD and Case 2 are 
compared in Figure 5.3. The DSMC shock standoff distance has been reduced as compared 
to Case 1 and now the location and magnitude of the peak translational temperatures of the 
CFD and DSMC Case 2 solutions compare quite well. The differences between the solutions 
are the larger equilibrium temperature within the shock layer and the increased electronic 
temperature that persists upstream of the location of maximum translational temperature, 
which is discussed below. 
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Figure 5.3: Comparison of temperature profiles between the CFD solution 
and Case 2 (electronic energy levels added). 


Comparisons of species mole fractions between the CFD solution and Case 2 are presented in 
Figure 5.4 and Figure 5.5 for neutral and charged species, respectively. In general, the 
species fractions compare well. As with the vibrational and rotational temperatures, there is 
noise in the mole fractions extracted from the DSMC simulation below around 0.05% 
because of the low sample size of these species. The DSMC solution contains trace species 
further upstream of the shock front as compared to the CFD solution. This is most likely 
caused by rarefaction effects where there are insufficient collisions to keep all particles 
within the shock layer. This also accounts for the increased electronic temperatures upstream 
of the maximum translational temperature as seen in Figure 5.3, and is analogous to what 
occurs with the vibrational temperature upstream of the shock, except with the electronic 
energy levels, not only do the molecules carry energy upstream of the shock, but also the 
atomic species. Even for a small population of excited particles, the electronic temperature 
can be large. 
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Figure 5.4: Comparison of neutral species mole fractions between the CFD 
solution and Case 2. 
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Figure 5.5: Comparison of charged species mole fractions between the CFD 
solution and Case 2. 

5.4.3 Quantum-Kinetic Solution 


As was discussed earlier, DSMC Case 3 switches chemistry modules from the traditional 
TCE model to the new Q-K model. The various temperature profiles from CFD and Case 3 
are compared in Figure 5.6. For this case, the peak translational temperature is slightly 
lower and moved upstream as compared to the CFD solution. This discrepancy may 
attributable both to the difference in the resultant equilibrium reaction rates discussed in 
Section 4 and also because the Q-K model does not rely on equilibrium-based information. 


68 


If the Q-K reaction rates are compared to the reaction rates used in the CFD solution (the first 
referenced rate in each figure) in Appendix C, there can be a large discrepancy, even though 
the Q-K rate is generally still well within the spread of reported reaction rates. 
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Figure 5.6: Comparison of temperature profiles between the CFD solution 
and Case 3 (with electronic energy levels and Q-K chemistry). 

Comparisons of species mole fractions are presented between the CFD solution and DSMC 
Case 3 in Figure 5.7 and Figure 5.8 for neutral and charged species, respectively. Again, 
the mole fractions generally compare well between the CFD and Case 2 solutions with some 
notable differences. One difference is again the presence of trace species well upstream of 
the shock layer. 
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Figure 5.7: Comparison of neutral species mole fractions between the CFD 
solution and Case 3. 
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Figure 5.8: Comparison of charged species mole fractions between the CFD 
solution and Case 3. 


5.4.4 Comparison of Surface Heating 


The FIRE II calorimeter data presented by Cornette[21] provide heating values that contain 
the convective heating plus the contribution of the radiative flux absorbed by the beryllium 
calorimeter. In order to compare to these data, the code HARA was used to post-process the 
flow field solutions to obtain the non-equilibrium radiation heating to the surface. It was 
detennined[21] that the fraction of the radiative flux absorbed by the calorimeter, a, was 
0 . 72 . 
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The computed heat flux to the surface of the vehicle is compared to the values inferred from 
the flight data at a radial location of approximately 0.02 m in Table 5.3. When comparing 
the convective heating component between the CFD and fully catalytic DSMC solutions 
including electronic energy levels, very good agreement is found. The solutions agree to 
within 1% of each other. However, DSMC Case 1, without electronic energy levels, has a 
convective heating rate approximately 15 W/cm 2 higher, again highlighting the fact that 
electronic energy levels are required for the simulation. When the cases with electronic 
energy levels are compared to the inferred flight value, the convective heating is over- 
predicted by 34%. If the partially catalytic solutions are then considered, they under-predict 
the convective heating by around 15%. Therefore, the surface of the flight vehicle was likely 
at least partially catalytic to atomic recombination. 

The radiative heating component is also listed in Table 5.3 in the third column. Only the 
DSMC simulations with electronic energy levels modeled have been post-processed for 
radiative heating. Because of the differences in molar fractions between the DSMC 
simulations and the CFD solution, there is about a 10 W/cm 2 increase in radiative heating for 
the DSMC simulations. There are generally more species present in the shock layer that 
radiate in the DSMC simulations, such as N, O, and N 2 + . The increase in radiative heating 
for the partially catalytic DSMC simulations is accounted for by the increased levels of N 
and O near the surface since these species do not recombine upon colliding with the wall. 


Table 5.3: FIRE II surface heating at r = 0.02 m. 


Case 

Convective Heating (W/cm 2 ) 

Radiative Heating (W/cm 2 ) 

Qconv CXC|rad 

(W/cm 2 ) 

CFD 

217 

25 

235 

Case 1 

230 

— 

— 

Case 1 pc 

148 

— 

— 

Case 2 

219 

36 

245 

Case 2 pc 

133 

95 

201 

Case 3 

217 

35 

242 

Case 3 pc 

137 

72 

189 

Flight 

162 

18 

175 
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Section 6 
Conclusions 

This section contains a summary of the results of this thesis, the original contributions made 
to the field, and identifies future research directions. 

6.1 Summary of Results 

Extensions to the recently introduced kinetic-theory approach for computing chemical 
reaction rates (the Q-K method) have been developed and validated for DSMC for electronic 
energy level transitions and reactions involving charged species. These new models make 
use of the principles of microscopic reversibility and molecular-level energy exchange to 
predict the probability that an energy level transition or a chemical reaction occurs during a 
collision between two particles. Procedures have been defined that are required to implement 
these models in the DSMC technique. 

For each of these new Q-K reaction-rate models, comparisons have been made between the 
predicted analytic curves, the sampled transition or reaction rates from DSMC simulations, 
and rates obtained from the literature. The analytic and sampled values compare favorably to 
the rates from the literature; the historical rates often differ by several orders of magnitude 
for the same reactions. It should be noted again, however, that the Q-K model is 
phenomenological in nature and there may be cases for which it leads to rates that are outside 
of the range of uncertainty of the measured rates. Even so, the usefulness of these models 
has been shown as engineering models, especially where there are no data or the uncertainty 
in reported values is quite large. If there are reactions for which there is a high degree of 
certainty in the reported rate, they can be specified to use the TCE model. 

The proposed extensions to the Q-K model were implemented in the DSMC Analysis Code 
(DAC) of NASA and tested on the Project Fire II flight experiment. The addition of 
electronic energy levels to the standard chemistry model resulted in better agreement 
between DSMC and CFD. Both the location and magnitude of the maximum translational 
temperature agree. The DSMC solution without electronic energy levels had predicted a 
higher peak translational temperature and shock thickness approximately 40% greater than 
the CFD solution. However, because of rarefaction effects, trace species with elevated 
electronic energies were found upstream of the shock layer, resulting in an elevated 
electronic temperature. Similar results were found when the chemistry model was changed 
from the TCE model to the Q-K model. The only major difference was a slightly thicker 
shock layer and lower peak translational temperature because of the differences in resultant 
reaction rates. Most likely, the Q-K model better represents the reaction rates because it is a 
phenomenological approach, whereas for the other models, the reaction rates are extrapolated 
from low temperature experiments. 

Once surface convective heating had been obtained, the flow field was post-processed using 
the HARA non-equilibrium radiation solver to compare DSMC, CFD and values obtained 
from the flight experiment. The convective heating of the fully catalytic DSMC simulations 
with electronic energy levels included compared to within 1% of values obtained using CFD. 
These values are, however, around 34% higher than the value inferred from flight. When a 
partially catalytic surface boundary condition is employed, allowing only electron/ion 
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recombinations, the convective heating dropped to levels on the order of 15% lower than the 
flight value. Therefore, the fully catalytic and partially catalytic simulations have been 
shown to bound the flight value of convective heating. The radiative heating computed from 
the DSMC simulations is around 10 W/cnr higher than that of the CFD simulation, which is 
caused by slightly larger amounts of N, O, and N 2 + . However, the nonequilibrium radiation 
solver is required to make certain assumptions in its calculations, so a direct input of the 
sampled non-Boltzmann distribution of electronic energy levels from the DSMC simulation 
would negate the necessity of the assumptions. 

6.2 Contributions 

The ability to simulate rarefied, ionized, hypersonic flows has been added to the suite of 
simulation capabilities for the National Aeronautics and Space Administration. Specifically, 
physical models for the treatment of charged particles using the DSMC method have been 
implemented in DAC. It has been shown that the inclusion of electronic energy levels, which 
few other DSMC solvers have the capability to do, is necessary for lunar-reentry velocity 
missions for the accurate prediction of shock shape and surface heating. The ability to 
predict electronic energy level transitions and resulting distributions using only energy level 
data has been developed as a result of this thesis. Also, the QK chemistry model was 
extended to include reactions involving charged species. This model allows the inclusion of 
chemical reactions for which specific reaction rate data is not known and may be more 
appropriate for non-equilibrium conditions. Both electronic energy level transition rates and 
chemical reaction rates have been shown to agree well with rates from the literature. 


73 



Section 7 
References 

1 . Bose, D., et al., Uncertainty Assessment of Hypersonic Aerothermodynamics 
Prediction Capability, in 42nd AIAA Thermophysics Conference 2011, AIAA Paper 
No. 2011-3141: Honolulu, HI. 

2. Bird, G.A., Molecular Dynamics and the Direct Simulation of Gas Flows. 1994, 
Oxford, UK: Oxford University Press. 

3. Ozawa, T., et al., Quasi-Classical Trajectory Modeling of OH Production in Direct 
Simulation Monte Carlo. Journal of Thennophysics and Heat Transfer, 2005. 19(2): 
p. 235-244. 

4. Farbar, E.D., Kinetic Simulation of Rarefied and Weakly Ionized Hypersonic Fow 
Fields, in Aerospace Engineering^) 1 0. The University of Michigan: Ann Arbor, 
MI. p. 170. 

5. Farbar, E.D. and I.D. Boyd, Simulation of Reactions Involving Charged Particles in 
Hypersonic Rarefied Flows, 2009, AIAA Paper No. 2009-267. 

6. Farbar, E.D. and I.D. Boyd, Simulation of FIREII II Reentry Flow Using the Direct 
Simulation Monte Carlo Method, 2008, AIAA Paper No. 2008-4103. 

7. Ozawa, T., et al., DSMC-CFD Comparison of a High Altitude, Extreme-Mach 
Number Reentry Flow, 2008, AIAA Paper No. 2008-1216. 

8. Ozawa, T., J. Zhong, and D.A. Levin, Development of kinetic-based energy 
exchange models for noncontinuum, ionized hypersonic flows. Physics of Fluids, 
2008. 20. 

9. Ozawa, T., et al., Modeling of the Stardust Reentry Flows with Ionization in DSMC, 
in 45th AIAA Aerospace Sciences Meeting and ExhibitlOOl , AIAA Paper No. 2007- 
611: Reno, NV. 

10. Bird, G.A. A Comparison of Collision Energy-based and Temperature-based 
Procedures in DSMC. in 26th Symposium on Rarefied Gas Dynamics. 2009. Kyoto, 
Japan: American Institue of Physics. 

11. Bird, G.A. Chemical Reactions in DSMC. in 27th Symposium on Rarefied Gas 
Dynamics. 2010. Pacific Grove, CA: American Institute of Physics. 

12. Bird, G.A. The Quantum-Kinetic Chemistry Model, in 27th International 
Symposium on Rarefied Gas Dynamics. 2010. Pacific Grove, CA. 

13. Bird, G.A., The Q-K Model for Gas-Phase Chemical Reaction Rates. Physics of 
Fluids, 2011. 23(10): p. 106101. 

14. Gallis, M.A., R.B. Bond, and J.R. Torczynski, A kinetic-theory approach for 
computing chemical-reaction rates in upper-atmosphere hypersonic flows. The 
Journal of Chemical Physics, 2009. 131(12): p. 124311. 

15. Bird, G.A., Nonequilibrium Radiation During Re-Entry at 10 km/s, 1987, AIAA 
Paper No. 1987-1543. 


74 



16. Gallis, M.A. and J.K. Harvey, Nonequilibrium Thermal Radiation from Air Shock 
Layers Modeled with Direct Simulation Monte Carlo. Journal of Thermophysics 
and Heat Transfer, 1994. 8(4): p. 765-772. 

17. Carlson, A.B. and H.A. Hassan, Radiation Modeling with Direct Simulation Monte 
Carlo. Journal of Thermophysics and Heat Transfer, 1992. 6: p. 631-636. 

18. Ozawa, T., et ah, Development of Coupled Particle Hypersonic Flowfield-Photon 
Monte Carlo Radiation Methods. Journal of Thermophysics and Heat Transfer, 
2010. 24(3): p.612-622. 

19. Johnston, C.O., B.R. Hollis, and K. Sutton, Spectrum Modeling for Air Shock-Layer 
Radiation at Lunar-Return Conditions. Journal of Spacecraft and Rockets, 2008. 
45(5): p. 865-878. 

20. Johnston, C.O., B.R. Hollis, and K. Sutton, Non-Boltzmann Modeling for Air 
Shock-Layer Radiation at Lunar-Return Conditions. Journal of Spacecraft and 
Rockets, '2008. 45(5): p. 879-890. 

2 1 . Cornette, E.S., Forebody Temperature and Calorimeter Heating Rates Measured 
During Project Fire II Reentry at 11.35 km/s, 1966, NASA TM X-1305. 

22. Birdsall, C.K. and A.B. Langdon, Plasma Physics via Computer Simulation. 1985, 
New York: McGraw-Hill. 

23. Hockney, R.W. and J.W. Eastwood, Computer Simulation Using Particles. 1988, 
Bristol, U.K.: Inst, of Phys. Publishing. 

24. Bird, G.A., Molecular Gas Dynamics. 1976, Oxford, U.K.: Clarendon Press. 

25. Taylor, J.C., A.B. Carlson, and H.A. Hassan, Monte Carlo Simulation of Radiating 
Re-Entrv Flow. Journal of Thermophysics and Heat Transfer, 1994. 8(3): p. 478- 
485. 

26. Bird, G.A., Low Density Aerothermodynamics, 1985, AIAA Paper No. 1985-0994. 

27. Bird, G.A., Direct Simulation of Typical AOTV Entry Flows, 1986, AIAA Paper 
No. 1986-1310. 

28. Gallis, M.A. and J.K. Harvey, Ionization Reactions and Electric Fields in Plane 
Hypersonic Shock Waves. Progess in Astronautics and Aeronautics, 1992. 160 : p. 
234-244. 

29. Carlson, A.B. and H.A. Hassan, Direct Simulation of Reentry Flows with 
Ionization. Journal of Thermophysics and Heat Transfer, 1992. 6(3): p. 400-404. 

30. Boyd, I.D., Monte Carlo Simulation of Nonequilibrium Flow in a Low-power 
Hydrogen Arcjet. Physcis of Fluids, 1997. 9(10): p. 4575-4584. 

3 1 . Nanbu, K., Probability Theory of Electron-Molecule, Ion-Molecule, Molecule- 
Molecule, and Coulomb Collisions for Particle Modeling of Materials processing 
Plasmas and Gases. IEEE Transactions of Plasma Science, 2000. 28(3): p. 971— 
990. 


75 



32. Nanbu, K., Stochasitc Theory of Motion and Collision of Charged Particle in a 
Uniform Electric Field. Journal of the Physical Society of Japan, 1994. 63: p. 979- 
983. 

33. Nanbu, K., Direct Simulation Scheme Derived from the Boltzmann Equation. -II. 
Multicomponent Gas Mixtures. Journal of the Physical Society of Japan, 1980. 49: 
p. 2042-2049. 

34. Nanbu, K., Direct Simulation Scheme derived from the Boltzmann Equation-1: 
Monocomponent Gases. Journal of the Physical Society of Japan, 1980. 49: p. 
2042-2049. 

35. Babovsky, H. and R. Illner, A Convergence Proof for Nanbu' s Simulation Method 
for the Full Boltzmann Equation. SIAM J. Num. Anal., 1989. 26: p. 45-65. 

36. Wilmoth, R.G., G.J. LeBeau, and A.B. Carlson, DSMC Grid Methodologies for 
Computing Low-Densitv, Hvpersonic Flows About Reusable Launch Vehicles, 

1996, AIAA 1996-1812. 

37. LeBeau, G.J. and F.E. Lumpkin III, Application Highlights of the DSMC Analysis 
Code (DAC) Software for Simulating Rarefied Flows. Computer Methods in 
Applied Mechanics and Engineering, 2001. 191(6-7): p. 595-609. 

38. LeBeau, G.J., K.A. Boyles, and F.E. Lumpkin III, Virtual Sub-Cells for the Direct 
Simulation Monte Carlo Method, 2003, AIAA 2003-1031. 

39. Gallis, M. and J.R. Torczynski, Effect of Collision-Partner Selection Schemes on 
the Accuracy and Efficiency of the Direct Simulation Monte Carlo Method. 
International Journal for Numerical Methods in Fluids, 2010. 

40. Burt, J.M., E. Josyula, and I.D. Boyd, Techniques for Reducing Collision 
Separation in Direct Simulation Monte Carlo Calculations, in 42nd AIAA 
Thermophysics Conference 2011, AIAA Paper No. 2011-3312: Honolulu, HI. 

41. Bird, G.A., Monte Carlo simulation in an engineering context. Progess in 
Astronautics and Aeronautics, 1981. 74(1). 

42. Borgnakke, C. and P.S. Larsen, Statistical collision model for Monte Carlo 
simulation of polyatomic gas mixtures. Journal of Computational Physics, 1975. 18. 

43. Bergemann, F. and I.D. Boyd, New Discrete Vibrational Energy Model for Direct 
Simulation Monte Carlo Method. Progess in Astronautics and Aeronautics, 1994. 
160: p. 174-183. 

44. Haas, B.L., Thermochemistry Models Applicable to a Vectorized Particle 
Simulation, in Department fo Aeronautics and Astronautics 1990, Stanford 
University. 

45. Lumpkin III, F.E., B.L. Haas, and I.D. Boyd, Resolution of Differences Between 
Collision Number Definitions in Particle and Continuum Simulations. Physics of 
Fluids A, 1991. 3: p. 2282. 


76 



46. Haas, B.L., et al.. Rates of Thermal Relaxation in Direct Simulation Monte Carlo 
Methods. Physics of Fluids, 1994. 6(6): p. 2191-2201. 

47. Jeans, J.H., The Dynamical Theory of Gases. 2 ed. 1916, London: Cambridge 
University Press. 

48. Landau, L. and E. Teller, Zur Theorie der Schalldispersion. Phys. Z. Sowjetunion, 
1936. 10: p. 14. 

49. Gimelshein, S.F., I.D. Boyd, and M.S. Ivanov, Modeling of Internal Energy 
Transfer in Plume Flows of Polyatomic Molecules by the DSMC Method, 1999, 
AIAA Paper No. 1999-738. 

50. Gimelshein, N.E., S.F. Gimelshein, and D.A. Levin, Vibrational Relaxation Rates 
in the Direct Simulation Monte Carlo Method. Physics of Fluids, 2002. 14(12): p. 
4452-4455. 

5 1 . Boyd, I.D., Relaxation of Discrete Rotational Energy Distributions Using a Monte 
Carlo Method. Physics of Fluids A, 1993. 5(9): p. 2278-2286. 

52. Parker, J.G., Rotational and Vibrational Relaxation in Diatomic Gases. Physics of 
Fluids, 1959. 2(4): p. 449-462. 

53. Millikan, R.C. and D.R. White, Systematics of Vibrational Relaxation. Journal of 
Chemical Physics, 1963. 39(12): p. 3209-32)3. 

54. Gallis, M.A. and J.K. Harvey, Atomic Species Radiation from Air Modeled with 
Direct Simulation Monte Carlo Method. Journal of Thermophysics and Heat 
Transfer, 1995. 9(3): p. 456-463. 

55. Bird, G.A. Simulation of Multi-Dimensional and Chemically Reacting Flows, in 
11th International Symposium on Rarefied Gas Dynamics. 1978. Cannes, France: 
Commissariat a Lenergie Atomique. 

56. Tonks, L. and I. Langmuir, General Theory of a Plasma Arc. Physical Review, 
1929. 34(1): p. 876-922. 

57. Bartel, T.J. and C.R. Justiz, DSMC Simulation of Ionized Rarefied Flows, in 24th 
AIAA Fluid Dynamics Conference 1993, AIAA Paper 1993-3095: Orlando, FL. 

58. Justiz, C.R. and C. Dalton, A Hybrid Flow Model for Charged and Neutral 
Particles Around Spacecraft in Low Earth Orbit, in 27th AIAA Thermophysics 
Conference 1992, AIAA Paper 1992-2935: Nashville, TN. 

59. Gallis, M., R. Prasad, and J.K. Harvey, The Effect of Plasmas on the Aerodynamic 
Performance of Vehicles, 1998, AIAA 1998-2666. 

60. Park, C., Nonequilibrium Hypersonic Aerothermodynamics. 1990, New York: 
Wiley. 

6 1 . Lord, R.G., Modeling Vibrational Energy Exchange of Diatomic Molecules Using 
the Morse Interatomic Potential. Physics of Fluids, 1998. 10(3): p. 4. 


77 



62. Koura, K. and H. Matsumoto, Variable soft sphere molecular model for inverse- 
power-law or Lennar d- Jones potential. Physics of Fluids A, 1991. 3. 

63. Park, C., R.L. Jaffe, and H. Partridge, Chemical-Kinetic Parameters of Hyperbolic 
Earth Entry. Journal of Thermophysics and Heat Transfer, 2001. 15(1): p. 76-90. 

64. Kang, S.W. and M.G. Dunn, Theoretical and Measured Electron-Density 
Distributions for the RAM Vehicle at High Altitudes, 1972, AIAA Paper No. 1972- 
689. 

65. Gimelshein, N.E. and D.A. Levin, Modeling of Glow Radiation in the Rarefied 
Flow About an Orbiting Spacecraft. Journal of Thermophysics and Heat Transfer, 
2000. 14(4): p.471-479. 

66. Johnston, H.S., Gas Phase Reaction Kinetics of Neutral Oxygen Species, U.S.D.o. 
Commerce, Editor 1968, National Bureau of Standards: Washington, D.C. 

67. Boyd, I.D. and T. Gokcen, Computation of Axisymmetric and Ionized Hypersonic 
Flows Using Particle and Continuum Methods. AIAA Journal, 1994. 32(9): p. 
1828-1835. 

68. Bourdon, A. and P. Vervisch, Study of a Low-Pressure Nitrogen Plasma Boundary 
Layer Over a Metallic Plate. Physics of Plasmas, 1997. 4: p. 4144-4157. 

69. Park, C., Chemical-Kinetic Problems of Future NASA Missions, 1991, AIAA Paper 
No. 1991-464. 

70. Losev, S.A., et ah, Thermochemical Nonequilibrium Kinetic Models in Strong 
Shock Waves on Air, 1994, AIAA Paper No. 1994-1990. 

71. Fujita, K., T. Yamada, and N. Ishii, Impacts of Ablation Gas Kinetics on Hyperbolic 
Entry Radiative Heating, in 44th AIAA Aerospace Sciences Meeting2006, AIAA 
Paper No. 2006-1185: Reno, NV. 

72. Gupta, R.N., Navier-Stokes and Viscous Shock-Layer Solutions for Radiating 
Hypersonic Flows, 1987, AIAA Paper No. 1987-1576. 

73. Phelps, A.V., Cross Sections and Swarm Coefficients for Nitrogen Ions and 
Neutrals in N2 and Argon Ions and Neutrals in Arfor Energies from 0.1 eV to 10 
keV. J. Phys. Chem. Ref. Data, 1991. 20: p. 557. 

74. Park, C., Review of Chemical-Kinetic Problems of Future NASA Missions, I: Earth 
Entries. Journal of Thermophysics and Heat Transfer, 1993. 7: p. 385-398. 

75. Park, C., J.T. Howe, and R.L. Jaffe, Review of Chemical-Kinetic Problems of 
Future NASA Missions, II: Mars Entries. Journal of Thermophysics and Heat 
Transfer, 1994. 8(1): p. 9-23. 

76. Gupta, R.N., et ah, A Review of Reaction Rates and Thermodynamic and Transport 
Properties for an 1 1 -Species Air Model for Chemical and Thermal Non equilibrium 
Calculations to 30,000 K, 1990, NASA RP-1232. 

77. Boyd, I.D., Modeling of Plasma Formation in Rarefied Hvpersonic Entry Flows, 
2007, AIAA Paper No.' 2007-206. 


78 



78. Carlson, A.B., Direct Simulation Monte Carlo with Ionization and Radiation, in 
Department of Mechanical and Aerospace Engineering 1990, North Carolina Sate 
University: Raleigh, p. 82. 

79. Tayal, S.S., Effective Collision Strengths for Electron Impact ofNI. Atomic Data 
and Nuclear Data Tables, 2000. 76: p. 191-212. 

80. Frost, R.M., et ah, Calculated cross sections and measured rate coefficients for 
electron-impact excitation of neutral and singly ionized nitrogen. Journal of 
Applied Physics, 1998. 84(6): p. 2989-3003. 

8 1 . Zatsarinny, O. and S.S. Tayal, Electron Collisional Excitation Rates for O I using 
the B-Spline R-Matrix Approach. The Astrophysical Journal Supplement Series, 
2003. 148: p. 575-582. 

82. Bhatia, A.K. and S.O. Kastner, The Neutral Oxygen Spectrum. I. Collisionally 
Excited Level Populations and Line Intensities under Optically Thin Conditions. 

The Astrophysical Journal Supplement Series, 1995. 96: p. 325-341. 

83. Chemyi, G.G. and S.A. Losev, Development of Thermal Protection Systems for 
Interplanetary Flight, R.I.o. Mechanics, Editor 1999, International Science and 
Technology Center: Moscow. 

84. Telulet, P., J.P. Sarrette, and A.M. Gomes, Calculation of Electron-Impact Inelastic 
Cross-Sections and Rate Coefficients for Diatomic Molecules. Application to Air 
Molecules. Journal of Qualitative Spectroscopy and Radiative Transfer, 1999. 62: p. 
549-569. 

85. Gorelov, V.A., et ah, Experimental and Numerical Study of Nonequilibrium 
Ultraviolet NO andN2+ Emission in a Shock Layer. Journal of Thermophysics and 
Heat Transfer, 1998. 12(2): p. 172-179. 

86. Johnston, C.O., Nonequilibruim Shock-Layer Radiative Heating for Earth and Titan 
Entry, in Aerospace Engin eeringl 0 0 6 , Virginia Polytechnic Institute and State 
University: Blacksburg, VA. p. 245. 

87. Gorelov, V.A., et al. Nonequilibrium Ionization and Radiation Behind Shock Wave 
in Martian Atmosphere, in Third European Symposium on Aerothermodynamics for 
Space Vehicles. 1998. Noordwijk, The Netherlands: European Space Agency. 

88. Park, C., Rate Parameters for Electronic Excitation of Diatomic Molecules II. 
Heavy Particle-Impact Processes, 2008, AIAA Paper No. 2008-1446. 

89. Flagan, R.C. and J.P. Appleton, Excitation Mechanisms of the Nitrogen First- 
Positive and First-Negative Radiation at High-Temperature. Journal of Chemical 
Physics, 1972. 56(3): p. 1163-1173. 

90. Nagy, O., Excitation Cross-Sections of N2+ Molecular Ion by Electron Impact and 
the Vibrational Energy Levels of the Three Target States. Chemical Physics, 2003. 
286(1): p. 109-114. 


79 



91. Guerra, V. and J. Loureiro, Electron and heavy particle kinetics in a low-pressure 
nitrogen glow discharge. Plasma Sources Science & Technology, 1997. 6: p. 361— 
372. 

92. Raynaud, E., et al., Huygens Aerothermal Environment: Radiative Heating, in 3rd 
International Planetary Probe Workshop2005, ESA: Anavyssos, Greece. 

93. Magin, T.E., et al., Nonequilibrium Radiation Modeling for Huygens Entry, in 3rd 
International Planetary Probe Works hop2005: Anavyssos, Greece. 

94. Cheatwood, F.M. and P. Gnoffo, User's Manual for the Langley 
Aerothermodynamic Upwind Relaxation Algorithm (LAURA), April 1996, NASA 
TM 4674. 

95. Gnoffo, P., R.N. Gupta, and J.L. Shinn. Conservation Equations and Physical 
Models for Hvpersonic Air Flows in Thermal and Chemical Nonequilibrium, in 
NASA TP 2867. Feb. 1989. 

96. Roe, P.L., Approximate Riernarm Solvers, Parameter Vectors, and Difference 
Schemes. Journal of Computational Physics, 1981. 43(2): p. 357-372. 

97. Yee, H.C., On Symmetric and Upwind TVD Schemes, in NASA TM 553251986. 

98. Mcbride, B.J., M.J. Zehe, and S. Gordon, NASA Glenn Coefficients for Calculating 
Thermodynamic Properties of Individual Species, Sept. 2002, NASA TP 2002- 
211556/ 

99. Sutton, K. and P. Gnoffo, Multi- Comp orient Diffusion with Application to 
Computational Aerothermodynamics, 1998, AIAA Paper 98-2575. 

100. Ralchenko, Y. and e. al. NIST Atomic Spectra Database, Version 3.1.0. July 2006; 
Available from: http://phvsics.nist.gov/PhysRefData/ASD/index.html . 

101. Team, T.O.P., The Opacity Project, Vol. 1. 1995, Bristol and Philedelphia: Institute 
of Physics Publishing. 

102. Cunto, W. and e. al, TOPbase at the CDS. Astronomy and Astrophysics, 1993. 275 : 
p. L5-L8. 

103. Soon, W.H. and J.A. Kune, Nitrogen Plasma Continuum Emission Associated with 
N-(3P) and N-(1D) Ions. Physical Review A, 1990. 41 : p. 4531-4533. 

104. Chauveau, S., et al., Radiative Transfer in LTE Air Plasmas for Temperatures up to 
15,000 K. Journal of Qualitative Spectroscopy and Radiative Transfer, 2003. 77: p. 
113-130. 

105. Chambers, L.H., Predicting Radiative Heat Transfer in Thermochemical 
Nonequilibrium Flow Fields, 1994, NASA TM-4564. 

106. Laux, C.O., Optical Diagnostics and Radiative Emission of Air Plasmas, 1993, 
High Temperature Gas Dynamics Lab, Mechanical Engineering Dept., Rept. T-288, 
Stanford University. 


80 


107. Whang, T.J., et al., Franck-Condon Factors of the Transitions ofN2. Journal of 

Quantitative Spectroscopy and Radiative Transfer, 1996. 55: p. 335-344. 

108. Stahel, D., M. Leoni, and K. Dresslar, Nonadiabatic Representations of the IZu and 
IFIu States of the N2 Molecule. Journal of Chemical Physics, Sept. 1983. 79: p. 
2541-2558.' 

109. Chauveau, S., et al., Contributions of Diatomic Molecular Electronic Systems to 
Heated Air Radiation. Journal of Quantitative Spectroscopy and Radiative Transfer, 
2002. 72: p. 503v530. 

1 10. Park, C., Radiation Enhancement by Nonequilibrium in Earth's Atmosphere. 

Journal of Spacecraft and Rockets, 1985. 22: p. 27-36. 

111. Hash, D., et al., FIRE II Calculations for Hypersonic Nonequilibrium 
Aerothermodynamics Code Verification: DPLR, LAURA and US3D, 2007, AIAA 
2007-605. 


81 



Appendices 


Appendix A: 

Species and Chemistry Data 


Table A.1 : Parameters used in the VHS molecular model and characteristic 

temperatures. 


Species 

T ref (K) 

dref( m) 

0) 

Orot( K) 

@vib (K) 

&diss (K) 

@ion (K) 

o 2 

288.0 

3.96E-10 

0.77 

2.080 

2256.0 

59500.0 

140066.5 

n 2 

288.0 

4.07E-10 

0.74 

2.874 

3371.0 

11355.0 

180798.3 

0 

288.0 

3.00E-10 

0.75 

- 

- 

- 

158053.5 

N 

288.0 

3.00E-10 

0.75 

- 

- 

- 

168613.6 

NO 

288.0 

4.00E-10 

0.79 

2.080 

2719.0 

75500.0 

107457.8 

O2 

288.0 

3.96E-10 

0.77 

2.080 

2256.0 

59500.0 

- 

n 2 + 

288.0 

4.07E-10 

0.74 

2.874 

3371.0 

11355.0 

- 

0 + 

288.0 

3.00E-10 

0.75 

- 

- 

- 

- 

N + 

288.0 

3.00E-10 

0.75 

- 

- 

- 

- 

NO + 

288.0 

4.00E-10 

0.79 

2.080 

2719.0 

75500.0 

- 

e 

288.0 

1.00E-10 

0.70 

- 

- 

- 

- 


Table A.2: Internal energy relaxation parameters. 


Species 

^ rZot 

z 'Hc,A” >x P £r*) 

c 1 

c 2 

o 2 

1 

56.5 

153.5 

n 2 

1 

9.1 

220 

NO 

1 

9.1 

220 

o 2 + 

1 

56.5 

153.5 

n 2 + 

1 

9.1 

220 

NO + 

1 

9.1 

220 


Table A.3: Parameters used in new vibrational relaxation model. 



o 2 

n 2 

NO 

Z re f 

17900.0 

52500.0 

123000.0 

T re f (K) 

2256.0 

3371.0 

2791.0 

e d (K) 

59500.0 

113500.0 

75500.0 


Table A.4: Electronic energy levels of 0 2 . 


Level 

9 

Eel (J) 

State 

1 

3 

0.0000E+01 


2 

2 

1.5728E-19 
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Level 

9 

Eel (J) 

State 

3 

1 

2.621 1 E-19 

9 

4 

1 

6.5663E-19 


5 

6 

6.8912E-19 


6 

3 

7.0306E-19 


7 

3 

9.9267E-19 



Table A.5: Electronic energy levels of N 2 . 


Level 

9 

Eel (3) 

State 

1 

1 

0.0000E+01 


2 

3 

9.9727E-19 


3 

6 

1.1843E-18 


4 

6 

1.1881E-18 


5 

3 

1.3165E-18 


6 

1 

1.3538E-18 


7 

2 

1.3763E-18 


8 

2 

1.4483E-18 


9 

5 

1.5415E-18 


10 

1 

1.6925E-18 


11 

6 

1.7242E-18 


12 

6 

1.7707E-18 


13 

10 

1.8474E-18 


14 

6 

1.9388E-18 


15 

6 

2.0778E-18 



Table A.6: Electronic energy levels of O. 


Level 

g 

E (J) 

Core 

Config. 

Term 

1 

9 

1.55E-21 

2s 2 

■'4- 

Q_ 

CM 

3 P 

2 

5 

3.15E-19 

2s 2 

Q_ 

CM 

'D 

3 

1 

6.71 E-19 

CM 

(/) 

CM 

2p 4 

'S 

4 

5 

1.47E-18 

2s 3 2p 2 

3s 

CJ1 

C/) 

0 

5 

3 

1.53E-18 

2s 3 2p 2 

3s 

CO 

W 

0 

6 

15 

1.72E-18 

2s 

3p 

5p 
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g 

E (J) 

Core 

Config. 

Term 

7 

9 

1.76E-18 

2s 3 2p 2 

3p 

3 P 

8 

5 

1.90E-18 

2s 3 2p 2 

4s 

5 S° 

9 

3 

1.91E-18 

2s 3 2p 2 

4s 

CO 

05 

O 

10 

25 

1.94E-18 

2s 3 2p 2 

3d 

5 d° 

11 

15 

1.94E-18 

2s 3 2p 2 

3d 

3 d° 

12 

15 

1.97E-18 

2s 3 2p 2 

4p 

Sp 

13 

9 

1.98E-18 

2s 3 2p 2 

4p 

3 P 

14 

15 

2.01E-18 

2s 3 2p 2 

3s 

CO 

D 

o 

15 

5 

2.03E-18 

2s 3 2p 2 

5s 

5 s° 

16 

3 

2.03E-18 

2s 3 2p 2 

5s 

3 s° 

17 

5 

2.04E-18 

2s 3 2p 2 

3s 

1 D° 

18 

25 

2.04E-18 

2s 3 2p 2 

4d 

5 d° 

19 

15 

2.04E-18 

2s 3 2p 2 

4d 

3 d° 

20 

35 

2.05E-18 

2s 3 2p 2 

4f 

5 f 

21 

21 

2.05E-18 

2s 3 2p 2 

4f 

3 f 

22 

15 

2.06E-18 

2s 3 2p 2 

5p 

5p 


9 

2.06E-18 

2s 3 2p 2 

5p 

3 P 

24 

25 

2.09E-18 

2s 3 2p 2 

5d 

5 d° 


15 

2.09E-18 

2s 3 2p 2 

5d 

3 d° 


35 

2.09E-18 

2s 3 2p 2 

5f 

5 f 

27 

21 

2.09E-18 

2s 3 2p 2 

5f 

3 f 

28 

288 

2.12E-18 

2s 3 2p 2 

n=6 

- 


392 

2.14E-18 

2s 3 2p 2 

n=7 

- 

30 

512 

2.15E-18 

2s 3 2p 2 

n=8 

- 


648 

2.15E-18 

2s 3 2p 2 

n=9 

- 


800 

2.16E-18 

2s 3 2p 2 

n=10 

- 


Table A.7: Electronic energy levels of N. 



g 


Core 

Config. 

Term 

1 

4 

0.0000E+01 

2s 2 

2p 3 

4 S° 

2 

10 

3.8195E-19 

2s 2 

CO 

Q_ 

CM 

2 d° 

3 

6 

5.7287E-19 

2s 2 

2p 3 

2po 

4 

12 

1.6554E-18 

2s 2 2p 2 

3s 

4 p 
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Level 

9 

Eel (J) 

Core 

Config. 

Term 

5 

6 

1.7122E-18 

2s 2 2p 2 

3s 

2 P 

6 

12 

1.7507E-18 

2s 

2p 4 

4 P 

7 

2 

1.8589E-18 

2s 2 2p 2 

3p 

2 S° 

8 

20 

1.8839E-18 

2s 2 2p 2 

3p 

4 D° 

9 

12 

1.8973E-18 

2s 2 2p 2 

3p 

4pO 

10 

4 

1.9219E-18 

2s 2 2p 2 

3p 

4 s° 

11 

10 

1.9235E-18 

2s 2 2p 2 

3p 

2 d° 

12 

6 

1.9426E-18 

2s 2 2p 2 

3p 

2po 

13 

10 

1.9798E-18 

2s 2 2p 2 

3s 

2 d 

14 

12 

2.0598E-18 

2s 2 2p 2 

4s 

4 P 

15 

6 

2.0698E-18 

2s 2 2p 2 

4s 

2 P 

16 

6 

2.0784E-18 

2s 2 2p 2 

3d 

2 P 

17 

28 

2.0802E-18 

2s 2 2p 2 

3d 

4 f 

18 

14 

2.0828E-18 

2s 2 2p 2 

3d 

2 f 

19 

12 

2.0828E-18 

2s 2 2p 2 

3d 

4 P 

20 

20 

2.0859E-18 

2s 2 2p 2 

3d 

4 d 

21 

10 

2.0884E-18 

2s 2 2p 2 

3d 

2 d 

22 

2 

2.1 151E-18 

2s 2 2p 2 

4p 

2 s° 

23 

20 

2.1220E-18 

2s 2 2p 2 

4p 

4 d° 

24 

12 

2.1258E-18 

2s 2 2p 2 

4p 

-P 

u 

o 

25 

10 

2.1300E-18 

2s 2 2p 2 

4p 

2 d° 

26 

4 

2.1343E-18 

2s 2 2p 2 

4p 

4 s° 

27 

6 

2.1377E-18 

2s 2 2p 2 

4p 

2po 

28 

90 

2.1912E-18 

2s 2 2p 2 

4d 2 P,4d 4 F,4d 4 D,4d 4 F,4d 4 P,4d 2 D 

29 

126 

2.1946E-18 

2s 2 2p 2 

4f 4 D,4f 4 F,4f 4 G,4f 2 D,4f 2 F,4f 2 G 

30 

450 

2.2368E-18 

2s 2 2p 2 

n=5 

31 

648 

2.2703E-18 

2s 2 2p 2 

n=6 

32 

822 

2.2864E-18 

2s 2 2p 2 

n=7 

33 

1152 

2.2968E-18 

2s 2 2p 2 

n=8 

34 

1458 

2.3040E-18 

2s 2 2p 2 

n=9 
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Level 

g 

Eel (j) 

Core 

Config. 

Term 

35 

1800 

2.3091E-18 

2s 2 2p 2 

n=10 


Table A.8: Electronic energy levels of NO. 


Level 

9 

Ee, (J) 

State 

1 

4 

0.0000E+01 

x 2 n 

r 

2 

8 

7.5485E-19 


3 

2 

8.7218E-19 


4 

4 

9.1 1 15E-19 


5 

4 

9.5349E-19 


6 

4 

9.7336E-19 


7 

4 

1.0343E-18 


8 

2 

1.0533E-18 


9 

4 

1.1979E-18 


10 

2 

1.2032E-18 


11 

4 

1 .2269E-18 


12 

4 

1.2401E-18 


13 

2 

1.2410E-18 


14 

2 

1.2485E-18 


15 

2 

1 .2515E-18 


16 

4 

1.2580E-18 



Table A.9: Electronic energy levels of 0 2 + . 


Level 

9 

Eel (J) 

State 

1 

4 

0.0000E+01 


2 

8 

6.5380E-19 


3 

4 

8.0594E-19 


4 

6 

8.0650E-19 


5 

4 

8.6013E-19 


6 

2 

9.2966E-19 


7 

4 

9.8329E-19 


8 

4 

1.0568E-18 
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Level 

9 

Eel (J) 

State 

9 

4 

1.2177E-18 

o 2 n 

u 

10 

4 

1.2276E-18 


11 

8 

1.301 IE-18 


12 

4 

1.311 IE-18 


13 

2 

1.3243E-18 


14 

2 

1.3786E-18 


15 

4 

1.4302E-18 



Table A.10: Electronic energy levels of N 2 + . 


Level 

9 

Eel (J) 

State 

1 

2 

0.0000E+01 


2 

4 

1.821 IE-19 


3 

2 

5.0578E-19 


4 

4 

5.0654E-19 


5 

8 

8.2636E-19 


6 

8 

9.1377E-19 


7 

4 

1.0492E-18 


8 

4 

1.0528E-18 


9 

4 

1.1323E-18 


10 

4 

1.1621E-18 


11 

8 

1.191 9E-1 8 


12 

8 

1.2316E-18 


13 

4 

1.2713E-18 


14 

4 

1.2733E-18 


15 

2 

1.2831E-18 


16 

2 

1.3309E-18 


17 

4 

1.4302E-18 



Table A.1 1 : Electronic energy levels of 0 + . 


Level 

9 

Eel (J) 

Core 

Config. 

Term 

1 

4 

0.0000E+01 

2s 2 

2p a 

4 S° 

2 

10 

5.3270E-19 

2s 2 

2p a 

2 d° 

3 

6 

8.0386E-19 

2s 2 

2p a 

2pO 
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Table A.12: Electronic energy levels of N + . 


Level 

9 

Ee,(J) 

Core 

Config. 

Term 

1 

9 

0.0000E+01 

2s 2 

2p 2 

3p 

2 

5 

3.0424E-19 

2s" 

2p" 

’D 

3 

1 

6.4929E-19 

2s" 

2p" 

n s 

4 

15 

1.8324E-18 

2s 

2p 3 

3 d° 

5 

9 

2.1696E-18 

2s 

2p 3 

3pO 

6 

5 

2.8642E-18 

2s 

2p 3 

1 D° 

7 

12 

2.9610E-18 

2s"2p 

3s 

3pO 


Table A.13: Electronic energy levels of NO + . 


Level 

fir 

Ee,( J) 

State 

1 

1 

0.0000E+01 


2 

3 

1.0367E-18 


3 

6 

1 .1771 E-18 


4 

6 

1.2293E-18 


5 

3 

1.3457E-18 


6 

1 

1.3814E-18 


7 

2 

1.4194E-18 


8 

2 

1.4595E-18 



Table A.14: Endothermic chemical reaction rate coefficients (m 3 /molecule/s) 
used in the TCE model (the first rate listed for each reaction is the 
rate used in the simulations). k f = a7*exp(-E a /fcT) 


Number 

Reaction 

a 

b 

c 

Ea(J) 

E rx (J) 

Reference 

1 

O 2 + O 2 — ► 0 + 0 + O 2 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[631 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



5.3930E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelsheinf651 

2 

O 2 + N 2 — ► 0 + 0 + N 2 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[631 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



1.1980E-11 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelsheinf651 

3 

O 2 + 0— »0 + 0 + 0 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 



1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa [9] 



1.4980E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelsheinf651 



4.5700E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Johnstonf661 

4 

o 2 + n^o + o + n 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 



1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang [641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



5.9930E-12 

- 1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

5 

O 2 + NO ^ 0 + 0 + NO 

3.321 IE-09 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 


88 







Number 

Reaction 

a 

b 

c 

Ea( J) 

E rx (J) 

Reference 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Birdfl 51 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

6 

O 2 + 02 + — ► 0 + 0 + 02 + 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[631 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Birdfl 51 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



5.3930E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

7 

o 2 + n 2 + ->■ 0 + 0 + n 2 + 

3.321 IE-09 

-1.50 

1.5 

8.1995E-19 

-8.1995E-19 

Park[631 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



1.1980E-11 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

8 

o 2 + o + ^ 0 + 0 + o + 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 



1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Birdfl 51 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



1.4980E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 



4.5700E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Johnston[661 

9 

0 2 + N + ^ 0 + 0 + N + 

1.6605E-08 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 



1.3750E-10 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

10 

0 2 + NO + 0 + 0 + NO + 

3.321 IE-09 

-1.50 

1.0 

8.1995E-19 

-8.1995E-19 

Park[631 



4.5800E-1 1 

-1.00 


8.1970E-19 

-8.1970E-19 

Bird[151 



5.9780E-12 

-1.00 


8.2150E-19 

-8.2150E-19 

Kang[641 



5.9930E-12 

-1.00 


8.1970E-19 

-8.1970E-19 

Gimelshein[651 

11 

N 2 + O 2 — ► N + N + O 2 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Birdfl 51 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang [64] 



4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[65] 

12 

N 2 + N 2 — ► N + N + N 2 

1.1620E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[15] 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



1.1600E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[9] 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 

13 

n 2 + o^n + n + o 

4.9980E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



1.8500E-08 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

14 

n 2 + n^n + n + n 

4.9980E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[631 



1.8500E-08 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



6.7830E-08 

-1.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd [671 



6.9000E-08 

-1.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

15 

N 2 + NO — > N + N + NO 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[15] 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[91 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

16 

N 2 + 0 2 + ^ N + N + O 2 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[91 


89 





Number 

Reaction 

a 

b 

c 

Ea( J) 

E rx (J) 

Reference 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

17 

n 2 + n 2 + — > N + n + n 2 + 

1.1620E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[151 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[64] 



1.1600E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[91 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 

18 

N 2 + 0 + — > N + N + 0 + 

4.9980E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



1.8500E-08 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

19 

N 2 + N + — > N + N + N + 

4.9980E-08 

-1.60 

1.0 

1.5630E-18 

-1.5630E-18 

Park[631 



1.8500E-08 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



6.7830E-08 

-1.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



7.1400E-08 

-1.50 


1.5630E-18 

-1.5630E-18 

Boyd [671 



6.9000E-08 

-1.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

20 

N 2 + NO + ^ N + N + NO + 

1.1620E-08 

-1.60 

0.5 

1.5630E-18 

-1.5630E-18 

Park[631 



6.1700E-09 

-1.60 


1.5610E-18 

-1.5610E-18 

Bird[1 51 



3.1550E-13 

-0.50 


1.5600E-18 

-1.5600E-18 

Kang[641 



4.9800E-08 

-1.60 


1.5630E-18 

-1.5630E-18 

Ozawa[91 



7.9700E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Boyd [671 



3.1870E-13 

-0.50 


1.5610E-18 

-1.5610E-18 

Gimelshein[651 

21 

N 2 + e — > N + N + e 

9.9670E-27 

2.60 

0.5 

1.5630E-18 

-1.5630E-18 

Bourdon[681 



4.9800E-06 

-1.60 


1.5610E-18 

-1.5630E-18 

Ozawa[91 



1.4900E-05 

-1.60 


1.5610E-18 

-1.5610E-18 

Park[691 



1.3450E-10 

-1.28 


1.5630E-18 

-1.5630E-18 

Losev[701 

22 

NO + O 2 — ► N + 0 + O 2 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[64] 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

23 

NO + N 2 — > N + 0 + N 2 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang [64] 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

24 

NO + O^N + O + O 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Ozawa[91 



1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

25 

NO + N^N + O + N 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang [64] 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[69] 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

26 

NO + NO — > N + 0 + NO 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fu jita[7 1 1 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang [64] 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

27 

NO + 0 2 + ^ N + 0 + 0 2 + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[1 51 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 


90 





Number 

Reaction 

a 

b 

c 

Ea( j) 

E rx (J) 

Reference 



6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimeishein[651 

28 

NO + N 2 + — > N + 0 + N 2 + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fu jita[7 1 1 



3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[151 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa [9] 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



6.5900E-10 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

29 

NO + 0 + — > N + 0 + 0 + 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fu jita[7 1 1 



7.6600E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[151 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Ozawa[9] 



1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

30 

NO + N + — > N + 0 + N + 

7.6600E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fu jita[7 1 1 



7.6600E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[151 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[91 



1.8300E-13 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

31 

NO + NO + — > N + 0 + NO + 

3.8300E-13 

-0.50 

1.0 

1.0420E-18 

-1.0420E-18 

Fujita[7 1 1 



3.8300E-13 

-0.50 


1.0430E-18 

-1.0430E-18 

Bird[151 



6.4760E-10 

-1.50 


1.0420E-18 

-1.0420E-18 

Kang[641 



3.3200E-09 

-1.50 


8.2140E-19 

-8.2140E-19 

Ozawa[9] 



8.3000E-15 

0.00 


1.0420E-18 

-1.0420E-18 

Park[691 



1.3180E-08 

-1.50 


1.0430E-18 

-1.0430E-18 

Gimelshein[651 

32 

N + e~ — > N + + e' + e‘ 

5.8100E-08 

-1.00 

0.0 

2.3220E-18 

-2.3220E-18 

Park[631 



4.1500E+03 

-3.82 


2.3300E-18 

-2.3300E-18 

Park[601 



9.0000E-1 1 

-1.00 


1.5000E-18 

-1.5000E-18 

Ozawa [71 



2.9890E-17 

0.60 


2.3300E-18 

-2.3300E-18 

Losev[70] 



5.8100E-08 

-1.00 


2.3300E-18 

-2.3300E-18 

Boyd [671 



1.0000E-14 

0.00 


2.3300E-18 

-2.3300E-18 

Bird[1 51 



1.8270E+02 

-3.14 


2.3330E-18 

-2.3300E-18 

Kang[641 



8.4340E-14 

0.00 


2.3300E-18 

-2.3300E-18 

Farbar[51 

33 

0 + e — > 0 + + e + e 

1.5900E-08 

-1.00 

0.0 

2.1880E-18 

-2.1880E-18 

Park[63] 



6.4800E+03 

-3.78 


2.1900E-18 

-2.1900E-18 

Park[601 



3.0000E-10 

-1.00 


1.4000E-18 

-1.4000E-18 

Ozawa [71 



8.6350E-18 

0.68 


2.1810E-18 

-2.1810E-18 

Losev[701 



1.5900E-08 

-1.00 


2.1810E-18 

-2.1810E-18 

Boyd [671 



3.0000E-12 

0.00 


2.1810E-18 

-2.1810E-18 

Birdfl 51 



5.9780E+01 

-2.91 


2.1810E-18 

-2.1810E-18 

Kang[641 



1.0540E-14 

0.00 


2.1810E-18 

-2.1810E-18 

Farbar[51 

34 

N 2 + 0^ NO + N 

9.9630E-17 

0.10 

0.0 

5.2470E-19 

-5.2470E-19 

Fu jita[7 1 1 



5.3000E-17 

0.10 


5.1750E-19 

-5.1750E-19 

Birdfl 51 



1.1620E-16 

0.00 


5.2460E-19 

-5.2460E-19 

Kang [64] 



9.4500E-18 

0.42 


5.9280E-19 

-5.9280E-19 

Ozawaf9] 



1.0600E-12 

-1.00 


5.1750E-19 

-5.1750E-19 

Park[691 



1.1200E-16 

0.00 


5.1750E-19 

-5.1750E-19 

Gimelsheinf651 

35 

NO + 0 — > 0 2 + N 

1.3950E-17 

0.00 

0.0 

2.6790E-19 

-2.6790E-19 

Park[631 



3.6000E-22 

1.29 


2.2380E-19 

-2.2380E-19 

Birdfl 51 



5.3140E-21 

1.00 


2.7200E-19 

-2.7200E-19 

Kang[641 



5.2790E-21 

1.00 


2.7190E-19 

-2.7190E-19 

Gimelsheinf651 

36 

N + N ^ N 2 + + e‘ 

7.3060E-23 

1.50 

0.0 

9.3190E-19 

-9.3190E-19 

Park[631 



3.3870E-17 

0.00 


9.3200E-19 

-9.3200E-19 

Park[60] 



2.9800E-20 

0.77 


9.3400E-19 

-9.3400E-19 

Birdfl 51 



2.3250E-17 

0.00 


9.3610E-19 

-9.361 0E-1 9 

Kang[641 

37 

N + 0 — NO + + e' 

8.8010E-18 

0.00 

0.0 

4.4040E-19 

-4.4040E-19 

Park[631 



8.7660E-18 

0.00 


4.4000E-19 

-4.4000E-19 

Park[601 
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Number 

Reaction 

a 

b 

c 

Ea( J) 

E rx (J) 

Reference 



1.4600E-21 

1.00 


4.4000E-19 

-4.4000E-19 

Park[691 



2.5500E-20 

0.37 


4.4220E-19 

-4.4220E-19 

Bird[151 



2.3250E-24 

1.50 


4.4040E-19 

-4.4040E-19 

Kang[641 



1.4990E-20 

0.50 


4.4730E-19 

-4.4730E-19 

Gupta[721 

38 

0 + 0 — > 0 2 + + e 

1.1790E-27 

2.70 

0.0 

1.1 130E-18 

-1.1 130E-18 

Fujita[7 1 1 



1.8300E-17 

0.00 


1.1 100E-18 

-1 .1 100E-18 

Park[601 



6.4200E-22 

0.49 


1.1200E-18 

-1.1200E-18 

Bird[151 



2.6570E-13 

-0.98 


1.1 160E-18 

-1 .1 160E-18 

Kang[641 

39 

N + + N 2 ^N 2 + + N 

1.6610E-18 

0.50 

1.0 

1.6840E-19 

-1.6840E-18 

Park[631 



1.6700E-17 

-0.18 


1.6700E-19 

-1.6700E-19 

Bird[151 



2.3250E-26 

2.10 


1.6700E-19 

-1.6700E-19 

Phelps[731 



3.3540E-19 

0.81 


1.7950E-19 

-1.7950E-19 

Kang[641 

40 

0 + + N 2 -► N 2 + + 0 

1.51 10E-18 

0.36 

1.0 

3.1480E-19 

-3.1480E-19 

Park[741 



1.4940E-18 

0.36 


3.1480E-19 

-3.1480E-19 

Losev[701 



1.0550E-16 

-0.21 


3.0600E-19 

-3.0600E-19 

Bird[151 



5.6460E-1 1 

-2.00 


3.1750E-19 

-3.1750E-19 

Kang[641 

41 

0 + + NO — > N + + 0 2 

2.3240E-25 

1.90 

1.0 

3.6730E-19 

-3.6730E-19 

Park[741 

42 

NO + + N -> N 2 + + 0 

1.1960E-16 

0.00 

1.0 

4.9010E-19 

-4.1090E-19 

Park[741 



2.8300E-17 

0.40 


4.9010E-19 

-4.901 0E-1 9 

Bird[151 

43 

NO + + N -► 0 + + N 2 

5.6460E-17 

-1.08 

1.0 

1.7670E-19 

-1.7670E-19 

Park[741 

44 

NO + + 0 -> N + + 0 2 

1.6610E-18 

0.50 

1.0 

1.0660E-18 

-1.0660E-18 

Park[741 



2.2250E-17 

0.31 


1.0670E-18 

-1.0670E-18 

Kang[641 

45 

NO + + 0 ^ O2 + N 

1.1960E-17 

0.29 

1.0 

6.7100E-19 

-6.7100E-19 

Fujita[7 1 1 

46 

NO + + 0 2 — O2 + NO 

3.9850E-17 

0.41 

1.5 

5.4010E-19 

-5.401 0E-1 9 

Park[741 



3.9850E-17 

0.00 


4.5010E-19 

-4.501 0E-1 9 

Losev[701 



1.7200E-14 

-0.17 


4.4730E-19 

-4.4730E-19 

Bird[151 



2.9890E-15 

0.17 


4.5560E-19 

-4.5560E-19 

Kang[641 

47 

0 2 + + N ^ N + + 0 2 

1.4450E-16 

0.14 

1.0 

3.9490E-19 

-3.9490E-19 

Park[741 

48 

02 + + N 2 — i ► N2 + + O 2 

1.6440E-17 

0.00 

1.5 

5.6190E-19 

-5.6190E-19 

Park[74] 

49 

0 2 + + 0 — 0 + + 0 2 

6.6420E-18 

-0.09 

1.0 

2.4850E-19 

-2.4850E-19 

Park[751 



1.8900E-16 

-0.52 


2.5900E-19 

-2.5900E-19 

Bird[151 



4.8490E-12 

-1.11 


3.8650E-19 

-3.8650E-19 

Kang [641 


Table A.15: Exothermic chemical reaction rate coefficients (m 3 /molecule/s) 
used in the TCE model. k r = a7*exp(-£ a //r7) 


Number 

Reaction 

a 

b 

c 

Ea( j) 

Erx( J) 

Reference 

50 

O + O + T — > O 2 + T 

6.3050E-44 

-0.50 

0.0 

0.0000E+01 

8.1995E-19 

Farbar[61 



8.3000E-45 

-0.50 


0.0000E+01 

8.1995E-19 

Gupta[761 

51 

n + n + t^n 2 + t 

5.691 0E-40 

-1.50 

0.0 

0.0000E+01 

1.5630E-18 

Farbar[61 



3.0060E-44 

-0.50 


0.0000E+01 

1.5630E-18 

Gupta[761 

52 

N + O + T^NO + T 

1.5830E-43 

-0.50 

0.0 

0.0000E+01 

1.0420E-18 

Farbar[61 



2.7850E-40 

-1.50 


0.0000E+01 

1.0420E-18 

Gupta[761 

53 

NO + N -> N 2 + O 

2.4900E-17 

0.10 

0.0 

0.0000E+01 

5.1750E-19 

Carlson[781 



2.0200E-17 

0.10 


0.0000E+01 

5.9280E-19 

Birdfl 5] 



4.0600E-12 

-1.36 


0.0000E+01 

5.1750E-19 

Boyd [671 



2.4900E-17 

0.00 


0.0000E+01 

5.1750E-19 

Gimelshein[651 

54 

0 2 +N — > NO + O 

5.2000E-22 

1.29 

0.0 

4.9680E-20 

2.2380E-19 

Birdfl 51 



4.1300E-21 

1.18 


0.0000E+01 

5.5300E-20 

Ozawa [91 



4.6000E-15 

-0.55 


0.0000E+01 

2.2380E-20 

Boyd [671 



1.5980E-18 

0.50 


4.9680E-20 

2.7190E-19 

Gimelshein[651 

55 

N 2 + + e‘ ->■ N + N 

4.4830E-12 

-0.50 

0.0 

0.0000E+01 

9.3400E-19 

Park[601 



8.8800E-10 

-1.23 


0.0000E+01 

9.3400E-19 

Birdfl 51 



1.0000E-08 

-1.43 


0.0000E+01 

9.3400E-19 

Ozawa[71 



7.2740E-12 

-0.65 


0.0000E+01 

9.3400E-19 

Boyd [771 
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Number 

Reaction 

a 

b 

c 

Ea( J) 

Erx(J) 

Reference 



1.5700E-17 

0.85 


0.0000E+01 

9.3400E-19 

Boyd [671 



2.491 0E-08 

-1.50 


0.0000E+01 

9.3400E-19 

Kang [641 

56 

NO + + e' — > N + 0 

1 .4940E-1 1 

-0.65 

0.0 

0.0000E+01 

4.4220E-19 

Park[601 



4.0300E-09 

-1.63 


0.0000E+01 

4.4220E-19 

Bird[1 51 



2.0000E-06 

-2.05 


0.0000E+01 

4.4220E-19 

Ozawa[7] 



1.3210E-09 

-1.19 


0.0000E+01 

4.4220E-19 

Boyd [771 



2.2000E-13 

-0.19 


0.0000E+01 

4.4220E-19 

Boyd [671 



1.3280E-08 

-1.50 


0.0000E+01 

4.4220E-19 

Kang [641 



2.9890E-1 1 

-1.00 


0.0000E+01 

4.4730E-19 

Gupta[721 

57 

C> 2 + + e — > 0 + 0 

2.4900E-12 

-0.50 

0.0 

0.0000E+01 

1.1200E-18 

Park[601 



3.8300E-09 

-1.51 


0.0000E+01 

1.1200E-18 

Bird[1 51 



6.0000E-08 

-1.60 


0.0000E+01 

1.1200E-18 

Ozawa[71 



1.4530E-04 

-2.41 


0.0000E+01 

1.1200E-18 

Boyd [771 



9.2200E-15 

0.29 


0.0000E+01 

1.1200E-18 

Boyd [671 



1.3280E-08 

-1.50 


0.0000E+01 

1 .1 160E-18 

Gupta[721 

58 

N 2 + + N — > N + + N2 

2.3700E-18 

-0.52 

1.0 

0.0000E+01 

1.6700E-19 

Birdfl 51 



2.3430E-14 

-0.61 


0.0000E+01 

1.6700E-19 

Boyd [771 



4.6500E-19 

0.50 


0.0000E+01 

1.6700E-19 

Kang [641 

59 

N 2 + + 0 0 + + n 2 

1.7700E-17 

-0.21 

1.0 

0.0000E+01 

3.1480E-19 

Birdfl 51 



1.9780E-18 

0.11 


0.0000E+01 

3.1480E-19 

Boyd [771 



4.1 180E-1 1 

-2.20 


0.0000E+01 

3.1480E-19 

Kang [641 

60 

N + + 0 2 — ► 0 + + NO 

2.4430E-26 

2.10 

1.0 

0.0000E+01 

2.1 120E-19 

Boyd [771 

61 

N 2 + + 0 ^ NO + + N 

4.1000E-18 

0.40 

1.0 

0.0000E+01 

4.9010E-19 

Birdfl 51 



1.7440E-18 

0.30 


0.0000E+01 

4.9010E-19 

Boyd [771 

62 

0 + + N 2 ^ NO + + N 

3.9700E-18 

-0.71 

1.0 

0.0000E+01 

1.7670E-19 

Boyd [771 

63 

N + + 0 2 ^ NO + + 0 

3.0400E-18 

0.29 

1.0 

0.0000E+01 

1.0660E-18 

Boyd [671 



1.6610E-16 

0.00 


0.0000E+01 

1.0670E-18 

Kang [641 

64 

0 2 + + N NO + + 0 

8.9180E-13 

-0.97 

1.0 

0.0000E+01 

6.7090E-19 

Boyd [771 

65 

0 2 + + NO — ► NO + + 0 2 

3.9900E-17 

0.41 

1.5 

0.0000E+01 

4.5010E-19 

Ozawa[9] 



6.1950E-16 

-0.05 


0.0000E+01 

4.5010E-19 

Boyd [771 



4.4700E-15 

-0.17 


0.0000E+01 

4.5010E-19 

Carlson[781 



2.9890E-17 

0.50 


0.0000E+01 

4.5010E-19 

Kang [641 

66 

N + + 0 2 — 0 2 + + N 



1.0 




67 

N 2 + + O2 — ► 02 + + N2 

4.5890E-18 

-0.04 

1.5 

0.0000E+01 

5.6190E-19 

Boyd [771 

68 

0 + + 0 2 -> 0 2 + + 0 

1.8900E-16 

-0.52 

1.0 

0.0000E+01 

2.5900E-19 

Birdfl 51 



4.9930E-18 

-0.00 


0.0000E+01 

2.5900E-19 

Boyd [771 



4.6500E-19 

0.50 


0.0000E+01 

2.5900E-19 

Kang [641 



4.9900E-14 

0.00 


0.0000E+01 

2.5900E-19 

Boyd [671 


Table A.16: Coefficients required for the collision volume in Q-K 
recombination reactions (V co n = a7*V re r). 


Number 

Reaction 

a 

b 

50 

O + O + O 2 — - > O 2 + O 2 

0.040 

-1.30 


O + O + N 2 — * O 2 + N 2 

0.070 

-1.20 


O + O + O — y O 2 + O 

0.080 

-1.20 


o + o + n-^o 2 + n 

0.090 

-1.20 


O + O + NO — ► 0 2 + NO 

0.065 

-1.20 


O + O + C>2 + — > O 2 + 02 + 

0.040 

-1.30 


0 + 0 + n 2 + -*• o 2 + n 2 + 

0.070 

-1.20 


0 + 0 + o + ^ o 2 + o + 

0.080 

-1.20 


O + O + N + -> 0 2 + N + 

0.090 

-1.20 


O + O + NO + -► 0 2 + NO + 

0.065 

-1.20 


O + O + e — > O 2 + e 

0.000 

0.00 
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Number 

Reaction 

a 

b 

51 

N + N + O 2 — > N 2 + O 2 

0.150 

-2.05 


N + N + N 2 -*• N 2 + N 2 

0.090 

-2.10 


n + n + o^n 2 + o 

0.160 

-2.00 


n + n + n->n 2 + n 

0.170 

-2.00 


N + N + NO — > N 2 + NO 

0.170 

-2.10 


N + N + 0 2 + N 2 + 0 2 + 

0.150 

-2.05 


N + N + N 2 + -*• N 2 + N 2 + 

0.090 

-2.10 


N + N + 0 + -► N 2 + 0 + 

0.160 

-2.00 


N + N + N + — > N 2 + N + 

0.170 

-2.00 


N + N + NO + — > N 2 + NO + 

0.170 

-2.10 


N + N + e- — > N 2 + e- 

0.000 

0.00 

52 

N + 0 + O 2 — > NO + O 2 

0.300 

-1.90 


N + 0 + N 2 ^ NO + N 2 

0.400 

-2.00 


N + 0 + 0 — > NO + 0 

0.300 

-1.75 


N + 0 + N — » NO + N 

0.300 

-1.75 


N + 0 + NO — ► NO + NO 

0.150 

-1.90 


N + 0 + 0 2 + — > NO + 0 2 + 

0.300 

-1.90 


N + 0 + N 2 + — *• NO + N 2 + 

0.400 

-2.00 


N + 0 + 0 + -> NO + 0 + 

0.300 

-1.75 


N + 0 + N + — > NO + N + 

0.300 

-1.75 


N + 0 + NO + — > NO + NO + 

0.150 

-1.90 


N + 0 + e" — > NO + e' 

0.000 

0.00 


Table A.17: Coefficients required for the collision volume in Q-K exchange 
reactions (E a — a(T/7" re f) |Eexch| j E a — E a — E eX ch if E exC h ^ 0.0). 


Number 

Reaction 

a 

b 

Eexch (J) 

34 

N 2 + O — ► NO + N 

0.150 

0.00 

-5.2465E-19 

35 

NO + O — > 0 2 + N 

0.100 

0.68 

-2.6785E-19 

39 

N + + N 2 — > N 2 + + N 

0.400 

0.10 

-1.6844E-19 

40 

0 + + N 2 -*• N 2 + + 0 

0.500 

0.20 

-3.1479E-19 

41 

0 + + NO — > N + + 0 2 

0.100 

0.50 

-3.6725E-19 

42 

NO + + N — ► N 2 + + O 

0.050 

0.00 

-4.9013E-19 

43 

NO + + N — ► 0 + + N 2 

0.050 

0.70 

-1.7672E-19 

44 

NO + + O — > N + + 0 2 

0.000 

0.00 

-1.0659E-19 

45 

NO + + O — >• 0 2 + + N 

0.000 

0.00 

-6.7100E-19 

46 

NO + + 0 2 -> 0 2 + + NO 

0.000 

0.00 

-4.5009E-19 

47 

0 2 + + N -»• N + + 0 2 

0.020 

0.10 

-3.9487E-19 

48 

0 2 + N 2 — * N 2 + 0 2 

0.050 

0.55 

-5.6192E-19 

49 

0 2 + + 0 -»• 0 + + 0 2 

0.700 

0.15 

-2.4852E-19 

53 

NO + N N 2 + O 

0.025 

0.90 

5.2465E-19 

54 

0 2 + N -> NO + O 

0.100 

0.48 

2.6785E-19 

58 

N 2 + + N N + + N 2 

0.400 

0.25 

1.6844E-19 

59 

n 2 + + O -> 0 + + n 2 

0.600 

0.35 

3.1479E-19 

60 

N + + 0 2 -»• 0 + + NO 

0.150 

0.40 

3.6725E-19 
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Number 

Reaction 

a 

b 

Eexch (J) 

61 

N 2 + + 0 -» NO + + N 

0.020 

0.80 

4.9013E-19 

62 

0 + + N 2 ^ NO + + N 

0.000 

0.00 

1.7672E-19 

63 

N + + 0 2 — > NO + + 0 

0.000 

0.00 

1.0659E-18 

64 

0 2 + + N NO + + 0 

0.000 

0.00 

6.7100E-19 

65 

0 2 + + NO — > NO + + 0 2 

0.070 

0.50 

5.4009E-19 

66 

N + + 0 2 — » 0 2 + + N 

0.000 

0.00 

3.9487E-19 

67 

N 2 + 0 2 — > 0 2 + n 2 

0.200 

0.30 

5.6192E-19 

68 

0 + + 0 2 -> 0 2 + + 0 

0.500 

0.10 

2.4852E-19 
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Appendix B: 

Quantum-Kinetic Electronic Energy Level Transition Rates 




Temperature (K) 


Temperature (K) 


Figure B.1 : Transitions for O 2 + 02- Figure B.3: Transitions for O 2 + N 2 . 
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Temperature (K) 


Temperature (K) 


Figure B.2: Transitions for Q 2 + O. 


Figure B.4: Transitions for Q 2 + N. 
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Temperature (K) 


Temperature (K) 


Figure B.5: Transitions for 0 2 + Figure B.8: Transitions for N 2 + 0 2 . 

NO. 




Temperature (K) Temperature (K) 


Figure B. 6: Transitions for N 2 + N 2 . Figure B. 9: Transitions for N 2 + O. 
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Figure B. 7: Transitions for N 2 + N. Figure B. 10: Transitions for N 2 + 

NO. 


97 



Temperature (K) 


Figure B.11: Transition from level 1- 
2 for N 2 + e. 



Figure B.14: Transition from level 2 
3 for N 2 + e. 



Temperature (K) 


Figure B.12: Transition from level 1- 
2 for O + e. 



Figure B.15: Transition from level 1 
3 for O + e. 



Temperature (K) 


Figure B.13: Transition from level 1- 
4 for O + e. 



Figure B.16: Transition from level 
1v5 for O + e. 
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Temperature (K) 


Figure B.17: Transition from level 2- 
3 for O + e. 



Figure B.20: Transition from level 2 
5 for O + e. 



Temperature (K) 


Figure B.18: Transition from level 1- 
2 for N + e. 



Figure B.21: Transition from level 1 
3 for N + e. 



Figure B.19: Transition from level 1- 
4 for N + e. 
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Temperature (K) 

Figure B.22: Transition from level 1 
5 for N + e. 
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Temperature (K) 


Figure B.23: Transition from level 2- 
3 for N + e. 
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Figure B.26: Transition from level 2 
4 for N + e. 
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Figure B.24: Transition from level 2- 
5 for N + e. 
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Figure B.27: Transition from level 3 
4 for N + e. 



Temperature (K) 


Figure B.25: Transition from level 3- 
5 for N + e. 



Figure B.28: Transitions for NO + 

o 2 . 
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Temperature (K) 



Figure B.29: Transitions for NO + Figure B.32: Transitions for NO + O. 

N 2 . 
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Figure B.30: Transitions for NO + N. Figure B.33: Transitions for NO + 

NO. 
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Figure B.31: Transition from level 1- 
3 for NO + e. 
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Figure B.34: Transition from level 1- 
4 for NO + e. 
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Temperature (K) 


Figure B.35: Transition from level 3- 
4 for NO + e. 



Temperature (K) 


Figure B.37: Transition from level 1- 
2 for N 2 + + e. 



Temperature (K) 


Figure B.36: Transition from level 1- 
3 for N 2 + + e. 



Temperature (K) 


Figure B.38: Transition from level 2- 
3 for N 2 + + e. 


Table B.1 : Sampled temperatures and their departure from the set 

equilibrium temperature. 


"^"equilibrium (K) 

10000.00 

% Difference 

20000.00 

% Difference 

30000.00 

% Difference 

Ttr (K) 

9949.81 

-0.502% 

19918.00 

-0.410% 

29901.50 

-0.328% 

"["electron (K) 

9953.79 

-0.462% 

19927.70 

-0.361% 

29942.00 

-0.193% 

Tel, 02 (K) 

10000.40 

0.004% 

20018.80 

0.094% 

30301.90 

1.006% 

T el, N2 (K) 

10011.80 

0.118% 

20017.00 

0.085% 

30060.00 

0.200% 

Tel,o(K) 

10019.80 

0.198% 

19992.40 

-0.038% 

29967.50 

-0.108% 

T e |, N (K) 

10009.50 

0.095% 

19996.20 

-0.019% 

30040.60 

0.135% 

T e |, NO (K) 

10011.20 

0.112% 

20037.70 

0.189% 

30119.40 

0.398% 
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T el, 02+ (K) 

10012.50 

0.125% 

20046.50 

0.233% 

30140.30 

0.468% 

T e |, N2+ (K) 

10015.70 

0.157% 

20031.50 

0.158% 

30106.20 

0.354% 

T e |, o+ (K) 

10011.30 

0.113% 

19993.80 

-0.031% 

29923.40 

-0.255% 

T el, N+ (K) 

10022.70 

0.227% 

19989.50 

-0.053% 

29896.20 

-0.346% 

Tel, NO+ (K) 

10011.20 

0.112% 

20014.50 

0.073% 

30062.80 

0.209% 


Table B.2: Sampled and Boltzmann distributions for 0 2 at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.7708E-01 

7.7719E-01 

-0.014% 

2 

1.6583E-01 

1.6585E-01 

-0.010% 

3 

3.8865E-02 

3.8807E-02 

0.149% 

4 

2.2448E-03 

2.2280E-03 

0.754% 

5 

1.0584E-02 

1.0566E-02 

0.175% 

6 

4.8032E-03 

4.7751 E-03 

0.588% 

7 

5.9326E-04 

5.8615E-04 

1.212% 


Table B.3: Sampled and Boltzmann distributions for 0 2 at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

5.5362E-01 

5.5315E-01 

0.086% 

2 

2.0892E-01 

2.0863E-01 

0.138% 

3 

7.1264E-02 

7.1363E-02 

-0.138% 

4 

1.7149E-02 

1.7099E-02 

0.292% 

5 

9.0780E-02 

9.1209E-02 

-0.470% 

6 

4.3126E-02 

4.3358E-02 

-0.535% 

7 

1.5139E-02 

1.5191 E-02 

-0.342% 


Table B.4: Sampled and Boltzmann distributions for 0 2 at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

4.2546E-01 

4.2474E-01 

0.170% 

2 

1.9476E-01 

1.9370E-01 

0.550% 

3 

7.5215E-02 

7.5191 E-02 

0.032% 

4 

2.9078E-02 

2.9007E-02 

0.245% 

5 

1.6010E-01 

1.6091 E-01 

-0.506% 

6 

7.7062E-02 

7.7793E-02 

-0.939% 

7 

3.8326E-02 

3.8662E-02 

-0.869% 
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Table B.5: Sampled and Boltzmann distributions for N 2 at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.9500E-01 

9.9505E-01 

-0.005% 

2 

2.1963E-03 

2.1776E-03 

0.859% 

3 

1.1353E-03 

1.1238E-03 

1.026% 

4 

1.1045E-03 

1.0937E-03 

0.992% 

5 

2.1789E-04 

2.1573E-04 

1.000% 

6 

5.5589E-05 

5.4867E-05 

1.316% 

7 

9.4670E-05 

9.3263E-05 

1.509% 


Table B.6: Sampled and Boltzmann distributions for N 2 at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.3201 E-01 

7.3248E-01 

-0.064% 

2 

5.9496E-02 

5.9350E-02 

0.246% 

3 

6.0508E-02 

6.0296E-02 

0.352% 

4 

5.9591 E-02 

5.9482E-02 

0.183% 

5 

1.8720E-02 

1.8681 E-02 

0.211% 

6 

5.4476E-03 

5.4391 E-03 

0.156% 

7 

1.0032E-02 

1.0029E-02 

0.034% 


Table B.7: Sampled and Boltzmann distribution for N 2 at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

3.6019E-01 

3.6023E-01 

-0.011% 

2 

9.7743E-02 

9.7283E-02 

0.473% 

3 

1.2460E-01 

1.2387E-01 

0.592% 

4 

1.2323E-01 

1.2275E-01 

0.391% 

5 

4.5125E-02 

4.5014E-02 

0.246% 

6 

1.371 IE-02 

1.371 IE-02 

0.001% 

7 

2.5906E-02 

2.5975E-02 

-0.265% 


Table B.8: Sampled and Boltzmann distributions for O at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.4529E-01 

9.4553E-01 

-0.025% 

2 

5.3797E-02 

5.3567E-02 

0.430% 

3 

8.1975E-04 

8.1264E-04 

0.875% 

4 

1.3074E-05 

1.2912E-05 

1.258% 
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5 

5.1522E-06 

5.01 19E-06 

2.799% 

6 

6.0933E-06 

6.0882E-06 

0.084% 

7 

2.7848E-06 

2.7387E-06 

1.682% 


Table B.9: Sampled and Boltzmann distributions for O at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.4519E-01 

7.4541 E-01 

-0.029% 

2 

1.3215E-01 

1.3224E-01 

-0.070% 

3 

7.2810E-03 

7.2842E-03 

-0.044% 

4 

2.0615E-03 

2.0531 E-03 

0.409% 

5 

9.9362E-04 

9.9083E-04 

0.282% 

6 

2.4553E-03 

2.4419E-03 

0.550% 

7 

1.2762E-03 

1.2686E-03 

0.597% 


Table B.10: Sampled and Boltzmann distributions for O at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

3.0872E-01 

3.0795E-01 

0.250% 

2 

8.0065E-02 

7.9930E-02 

0.170% 

3 

6.7504E-03 

6.7670E-03 

-0.245% 

4 

4.9942E-03 

4.9743E-03 

0.399% 

5 

2.5918E-03 

2.5813E-03 

0.406% 

6 

8.1065E-03 

8.0535E-03 

0.658% 

7 

4.4226E-03 

4.3898E-03 

0.748% 


Table B.11: Sampled and Boltzmann distributions for N at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

8.4636E-01 

8.4671 E-01 

-0.042% 

2 

1.3341 E-01 

1.331 IE-01 

0.225% 

3 

2.0090E-02 

2.0036E-02 

0.271% 

4 

1.5808E-05 

1.5762E-05 

0.292% 

5 

5.2667E-06 

5.2246E-06 

0.807% 

6 

7.9995E-06 

7.9046E-06 

1.201% 

7 

6.2563E-07 

6.0151 E-07 

4.011% 
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Table B.12: Sampled and Boltzmann distributions for N at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

4.3952E-01 

4.3890E-01 

0.141% 

2 

2.7547E-01 

2.7515E-01 

0.116% 

3 

8.2188E-02 

8.2688E-02 

-0.605% 

4 

3.2850E-03 

3.2799E-03 

0.155% 

5 

1.3379E-03 

1.3353E-03 

0.198% 

6 

2.3279E-03 

2.3227E-03 

0.223% 

7 

2.6346E-04 

2.6158E-04 

0.719% 


Table B.13: Sampled and Boltzmann distributions for N at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

1.0682E-01 

1.0620E-01 

0.580% 

2 

1.0632E-01 

1.0558E-01 

0.697% 

3 

3.9934E-02 

3.9954E-02 

-0.051% 

4 

5.8929E-03 

5.8551 E-03 

0.646% 

5 

2.5708E-03 

2.5524E-03 

0.722% 

6 

4.6881 E-03 

4.6517E-03 

0.783% 

7 

6.0134E-04 

5.9704E-04 

0.720% 


Table B.14: Sampled and Boltzmann distributions for NO at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.8591 E-01 

9.8601 E-01 

-0.010% 

2 

8.3771 E-03 

8.3267E-03 

0.605% 

3 

8.9642E-04 

8.8993E-04 

0.730% 

4 

1.3534E-03 

1.3421 E-03 

0.843% 

5 

9.9780E-04 

9.8763E-04 

1.030% 

6 

8.6417E-04 

8.5528E-04 

1 .040% 

7 

5.5615E-04 

5.5025E-04 

1.073% 


Table B.15: Sampled and Boltzmann distributions for NO at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.3783E-01 

7.3899E-01 

-0.157% 

2 

9.6384E-02 

9.6040E-02 

0.359% 

3 

1.5806E-02 

1.5699E-02 

0.684% 

4 

2.7450E-02 

2.7264E-02 

0.683% 

5 

2.3546E-02 

2.3388E-02 

0.675% 
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6 

2.1886E-02 

2.1765E-02 

0.557% 

7 

1.7534E-02 

1.7457E-02 

0.439% 


Table B.16: Sampled and Boltzmann distributions for NO at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

4.7041 E-01 

4.7202E-01 

-0.342% 

2 

1.5317E-01 

1.5259E-01 

0.383% 

3 

2.9009E-02 

2.8737E-02 

0.948% 

4 

5.2855E-02 

5.231 IE-02 

1.040% 

5 

4.7618E-02 

4.7228E-02 

0.826% 

6 

4.5222E-02 

4.5016E-02 

0.457% 

7 

3.8913E-02 

3.8862E-02 

0.132% 


Table B.17: Sampled and Boltzmann distributions for 0 2 + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.7135E-01 

9.7155E-01 

-0.020% 

2 

1.7156E-02 

1.7058E-02 

0.577% 

3 

2.8567E-03 

2.8334E-03 

0.821% 

4 

4.2667E-03 

4.2331 E-03 

0.795% 

5 

1.9302E-03 

1.9136E-03 

0.867% 

6 

5.8441 E-04 

5.7826E-04 

1.063% 

7 

7.9150E-04 

7.8423E-04 

0.927% 


Table B.18: Sampled and Boltzmann distributions for 0 2 + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

6.6647E-01 

6.6768E-01 

-0.181% 

2 

1.2558E-01 

1.2512E-01 

0.372% 

3 

3.6258E-02 

3.6057E-02 

0.557% 

4 

5.4176E-02 

5.3977E-02 

0.369% 

5 

2.9739E-02 

2.9632E-02 

0.361% 

6 

1.1570E-02 

1.1518E-02 

0.451% 

7 

1.9025E-02 

1.8970E-02 

0.293% 


Table B.19: Sampled and Boltzmann distributions for 0 2 + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

4.1 158E-01 

4.1287E-01 

-0.313% 

2 

1.7107E-01 

1.7034E-01 

0.427% 

3 

5.9420E-02 

5.8989E-02 

0.731% 
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4 

8.8882E-02 

8.8364E-02 

0.586% 

5 

5.1954E-02 

5.1754E-02 

0.386% 

6 

2.1926E-02 

2.1879E-02 

0.217% 

7 

3.8413E-02 

3.8443E-02 

-0.077% 


Table B.20: Sampled and Boltzmann distributions for N 2 + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

6.1227E-01 

6.1286E-01 

-0.096% 

2 

3.2813E-01 

3.2776E-01 

0.113% 

3 

1.5790E-02 

1.5718E-02 

0.457% 

4 

3.1316E-02 

3.1262E-02 

0.173% 

5 

6.2032E-03 

6.1665E-03 

0.596% 

6 

3.2954E-03 

3.2741 E-03 

0.649% 

7 

6.1995E-04 

6.1403E-04 

0.964% 


Table B.21: Sampled and Boltzmann distributions for N 2 + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

3.1457E-01 

3.1438E-01 

0.059% 

2 

3.2567E-01 

3.2514E-01 

0.163% 

3 

5.0346E-02 

5.0348E-02 

-0.003% 

4 

1.0017E-01 

1.0042E-01 

-0.245% 

5 

6.3131 E-02 

6.3071 E-02 

0.096% 

6 

4.5969E-02 

4.5958E-02 

0.025% 

7 

1.4075E-02 

1.4073E-02 

0.014% 


Table B.22: Sampled and Boltzmann distributions for N 2 + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

1.8259E-01 

1.8224E-01 

0.195% 

2 

2.3563E-01 

2.3481 E-01 

0.349% 

3 

5.3865E-02 

5.3742E-02 

0.229% 

4 

1.0748E-01 

1.0728E-01 

0.183% 

5 

9.9579E-02 

9.9136E-02 

0.447% 

6 

8.0402E-02 

8.0276E-02 

0.157% 

7 

2.8918E-02 

2.8946E-02 

-0.098% 
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Table B.23: Sampled and Boltzmann distributions for 0 + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.4566E-01 

9.4590E-01 

-0.025% 

2 

5.01 10E-02 

4.9904E-02 

0.414% 

3 

4.2276E-03 

4.2007E-03 

0.641% 


Table B.24: Sampled and Boltzmann distributions for 0 + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

6.9240E-01 

6.9214E-01 

0.037% 

2 

2.5131 E-01 

2.5137E-01 

-0.023% 

3 

5.6287E-02 

5.6491 E-02 

-0.361% 


Table B.25: Sampled and Boltzmann distributions for 0 + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

5.2579E-01 

5.2459E-01 

0.229% 

2 

3.6206E-01 

3.6242E-01 

-0.099% 

3 

1.121 6E-01 

1.1299E-01 

-0.735% 


Table B.26: Sampled and Boltzmann distributions for N + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.4103E-01 

9.4131 E-01 

-0.030% 

2 

5.8007E-02 

5.7735E-02 

0.471% 

3 

9.5709E-04 

9.4864E-04 

0.891% 

4 

2.7276E-06 

2.7013E-06 

0.974% 

5 

1.3391E-07 

1.4098E-07 

-5.016% 

6 

0.0000E+01 

5.1 142E-10 

-100.000% 

7 

9.9063E-10 

6.0919E-10 

62.615% 


Table B.27: Sampled and Boltzmann distributions for N + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

8.3493E-01 

8.3487E-01 

0.008% 

2 

1.5404E-01 

1.541 IE-01 

-0.046% 

3 

8.8414E-03 

8.8344E-03 

0.079% 
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4 

1.8248E-03 

1.8258E-03 

-0.056% 

5 

3.2492E-04 

3.2310E-04 

0.565% 

6 

1.4622E-05 

1.4505E-05 

0.810% 

7 

2.4824E-05 

2.4524E-05 

1.222% 


Table B.28: Sampled and Boltzmann distributions for N + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.6006E-01 

7.5955E-01 

0.067% 

2 

2.0205E-01 

2.0243E-01 

-0.188% 

3 

1.7556E-02 

1.7600E-02 

-0.251% 

4 

1.5094E-02 

1.5173E-02 

-0.519% 

5 

4.0278E-03 

4.0337E-03 

-0.145% 

6 

4.1872E-04 

4.1885E-04 

-0.031% 

7 

7.9535E-04 

7.9589E-04 

-0.068% 


Table B.29: Sampled and Boltzmann distributions for NO + at 10,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

9.9599E-01 

9.9603E-01 

-0.004% 

2 

1.6518E-03 

1.6379E-03 

0.849% 

3 

1.1971 E-03 

1.1854E-03 

0.988% 

4 

8.2190E-04 

8.1217E-04 

1.198% 

5 

1.7649E-04 

1.7472E-04 

1.015% 

6 

4.5549E-05 

4.4971 E-05 

1.285% 

7 

6.9164E-05 

6.8332E-05 

1.218% 


Table B.30: Sampled and Boltzmann distributions for NO + at 20,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

7.8349E-01 

7.8358E-01 

-0.012% 

2 

5.5182E-02 

5.5037E-02 

0.263% 

3 

6.6323E-02 

6.6215E-02 

0.163% 

4 

5.4691 E-02 

5.4809E-02 

-0.215% 

5 

1.7924E-02 

1.7975E-02 

-0.286% 

6 

5.2648E-03 

5.2652E-03 

-0.008% 

7 

9.1870E-03 

9.1786E-03 

0.091% 
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Table B.31: Sampled and Boltzmann distributions for NO + at 30,000 K. 


Level 

Sampled 

Boltzmann 

% Difference 

1 

4.5955E-01 

4.5873E-01 

0.178% 

2 

1.1342E-01 

1.1263E-01 

0.703% 

3 

1.6075E-01 

1.6052E-01 

0.145% 

4 

1.4048E-01 

1.4151E-01 

-0.728% 

5 

5.2954E-02 

5.3415E-02 

-0.862% 

6 

1.6242E-02 

1.6335E-02 

-0.568% 

7 

2.9677E-02 

2.9810E-02 

-0.446% 
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Appendix C: 

Quantum-Kinetic Reaction Rates 



Temperature (K) 


Figure C.1 : Q-K reaction rates for 
reaction 1. 



Temperature (K) 


Figure C.3: Q-K reaction rates for 
reaction 2. 



Temperature (K) 

Figure C.2: Q-K reaction rates for 
reaction 3. 



Temperature (K) 

Figure C.4: Q-K reaction rates for 
reaction 4. 
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Temperature (K) 

Figure C.5: Q-K reaction rates for 
reaction 5. 
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Figure C.8: Q-K reaction rates for 
reaction 6. 



Temperature (K) 


Figure C.6: Q-K reaction rates for 
reaction 7. 
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Figure C.9: Q-K reaction rates for 
reaction 8. 
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Figure C. 7: Q-K reaction rates for 
reaction 9. 



Temperature (K) 


Figure C.10: Q-K reaction rates for 
reaction 10. 
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Temperature (K) 

Figure C.11: Q-K reaction rates for 
reaction 11. 
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Figure C.12: Q-K reaction rates for 
reaction 13. 
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Figure C.13: Q-K reaction rates for 
reaction 15. 
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Figure C.14: Q-K reaction rates for 
reaction 12. 
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Figure C.15: Q-K reaction rates for 
reaction 14. 
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Figure C.16: Q-K reaction rates for 
reaction 16. 
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0 5xl0 3 lxio 4 2xl0 4 2xl0 4 3xl0 4 3xl0 4 4xl0 4 

Temperature (K) 

Figure C.17: Q-K reaction rates for 
reaction 17. 
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Figure C.20: Q-K reaction rates for 
reaction 18. 
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Figure C.18: Q-K reaction rates for 
reaction 19. 
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Figure C.21: Q-K reaction rates for 
reaction 20. 
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Figure C.19: Q-K reaction rates for 
reaction 21. 
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Figure C.22: Q-K reaction rates for 
reaction 22. 
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Figure C.23: Q-K reaction rates for 
reaction 23. 
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Figure C.26: Q-K reaction rates for 
reaction 24. 
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Figure C.24: Q-K reaction rates for 
reaction 25. 
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Figure C.27: Q-K reaction rates for 
reaction 26. 
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Figure C.25: Q-K reaction rates for 
reaction 27. 
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Figure C.28: Q-K reaction rates for 
reaction 28. 
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Figure C.29: Q-K reaction rates for 
reaction 29. 
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Figure C.32: Q-K reaction rates for 
reaction 30. 
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Figure C.30: Q-K reaction rates for 
reaction 31. 
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Figure C.33: Q-K reaction rates for 
reaction 32. 
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Figure C.31: Q-K reaction rates for 
reaction 33. 
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Figure C.34: Q-K reaction rates for 
reaction 34. 
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Figure C.35: Q-K reaction rates for 
reaction 35. 
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Figure C.38: Q-K reaction rates for 
reaction 36. 
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Figure C.36: Q-K reaction rates for 
reaction 37. 



Figure C.37: Q-K reaction rates for 
reaction 39. 
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Figure C.39: Q-K reaction rates for 
reaction 38. 
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Figure C.40: Q-K reaction rates for 
reaction 40. 
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Figure C.41: Q-K reaction rates for 
reaction 41 . 



Figure C.44: Q-K reaction rates for 
reaction 42. 



Figure C.42: Q-K reaction rates for 
reaction 43. 
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Figure C.45: Q-K reaction rates for 
reaction 44. 



Figure C.43: Q-K reaction rates for 
reaction 45. 
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Figure C.46: Q-K reaction rates for 
reaction 46. 
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Figure C.47: Q-K reaction rates for 
reaction 47. 



Figure C.50: Q-K reaction rates for 
reaction 48. 



Figure C.48: Q-K reaction rates for 
reaction 49. 
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Figure C.51: Q-K reaction rates for 
reaction 50. 
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Figure C.49: Q-K reaction rates for 
reaction 51. 
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Figure C.52: Q-K reaction rates for 
reaction 52. 
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Figure C.53: Q-K reaction rates for 
reaction 53. 
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Figure C.56: Q-K reaction rates for 
reaction 54. 
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Figure C.54: Q-K reaction rates for 
reaction 55. 
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Figure C.55: Q-K reaction rates for 
reaction 57. 
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Figure C.57: Q-K reaction rates for 
reaction 56. 
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Figure C.58: Q-K reaction rates for 
reaction 58. 
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Figure C.59: Q-K reaction rates for 
reaction 59. 



Figure C.62: Q-K reaction rates for 
reaction 60. 



Figure C.60: Q-K reaction rates for 
reaction 61. 
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Figure C.63: Q-K reaction rates for 
reaction 62. 
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Figure C.61: Q-K reaction rates for 
reaction 63. 
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Figure C.64: Q-K reaction rates for 
reaction 64. 
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Figure C.65: Q-K reaction rates for 
reaction 65. 



Figure C.67: Q-K reaction rates for 
reaction 66. 



Figure C.66: Q-K reaction rates for 
reaction 67. 
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Figure C.68: Q-K reaction rates for 
reaction 68. 
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